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is  compared  to  coda  from  deep  tel eseisms,  recorded  at  State  College,  Pennsylvania,  and 
it  is  seen  that  scattering  is  more  severe  at  PAS,  as  reflected  in  higher  coda  levels 
and  slower  decay  rate.  Consideration  of  energy  partitioning  and  coda  amplitude  sugges 
that  much  of  the  coda  consists  of  scattered  surface  waves.  Analysis  of  a  major 
Ps  conversion  arriving  3  s  after  direct  P  indicates  that  a  major  crustal  discontinuity 
at  about  20  km  depth  dips  at  moderate  angles  to  the  north  under  the  San  Gabriel 
Mountains.  This  interface  probably  represents  the  crustal  tectonic  boundary  between 
the  Transverse  Ranges  and  the  Los  Angeles  Basis. 

The  codas  of  long-period  Rayleigh  waves  recorded  at  WWSSN  and  Canadian  network 
stations  in  Western  North  America  from  eight  underground  explosions  at  NTS  are 
examined  in  an  effort  to  separate  scattering  and  anelastic  attenuation  effects.  Coda 
behavior  of  0.1  and  0.2  hz  Rayleigh  waves  follows  coda  characteristics  seen  in 
studies  of  short-period  S  waves.  Coda  decay  rate  is  seen  to  be  a  stable  observation 
over  most  stations  in  Western  North  America  and  is  consistent  with  the  hypothesis 
that  backscattered  surface  waves  from  heterogeneities  contained  within  the  western 
half  of  the  continent  form  the  Rayleigh  wave  coda.  The  basic  data  observables  of 
coda  level  and  decay  are  l'nfprnrpt ed  using  several  plausible  m_2clc.  The  single 
scattering  model  yields  a  coda  Q  consistent  with  previously  determined  Rayleigh 
anelastic  attenuation  coef f icients .  Separation  of  anelastic  and  scattering  Q  is 
possible  using  an  energy  flux  model  and  shows  that  scattering  Q  is  one  to  two  orders 
of  magnitude  higher  than  anelastic  Q.  However,  an  energy  flux  model  which  incorporate 
a  layer  of  scatterers  over  a  homogeneous  half-space  shows  that  all  Rayleigh  wave 
attenuation  can  be  explained  purely  by  scattering  effects  which  include  Rayleigh 
to  body  wave  conversion.  Coda  can  be  fit  equally  well  by  these  mutually  incompatible 
models.  It  is  not  likely  that  the  mechanisms  of  scattering  or  anelastic  attenuation 
can  be  addressed  by  coda  observations  of  a  single  homogeneous  data  set. 

P^^  waveforms  from  two  moderate-sized  earthquakes  in  Zambia  are  used  to  determine 
an  upper  mantle  P-wave  velocity  model  for  southern  Africa.  The  events  are: 

5/15/68,  m^=5.7,  depth  =  28  km;  and  12/2/68,  m^=5.9,  depth  *  6  km.  Focal  parameters 
for  these  events  are  constrained  by  previous  workers  from  teleseismic  body  wave 
inversion.  Synthetic  seismograms  are  generated  for  various  mantle  velocity  models 
using  a  wavenumber  integration  method  until  an  acceptable  fit  to  the  data  is  obtained. 
Quality  of  fit  is  measured  primarily  by  the  Pn/PL  amplitude  ratio.  Source-station 
geometry  also  allows  for  the  independent  sampling  of  the  upper  mantle  beneath  the 
Kapvaal-Rhodesian  craton  and  the  mobile  belt  provinces.  Synthetics  from  a  three-layer 
crust  over  half-space  mantle  model  do  not  show  prominent  precursor  arrivals 
seen  in  the  data;  these  are  interpreted  to  be  P-waves  turning  in  the  upper  mantle. 

The  synthetics  also  give  a  too  low  P^/PL  amplitude  ratio.  Synthetics  for  models 
with  a  mantle  P-wave  velocity  g  Jn-nt  of  0.00333/sec  fit  the  cratonic  path  data 
very  well,  since  there  is  no  in.  ion  of  interaction  with  a  low-velocity  zone, 
this  gives  a  minimum  lithospheric  i.ickness  of  120  km.  A  slightly  lower  gradient 
is  indicated  for  the  mobile  belt  regions,  with  a  minimum  lithospheric  thickness  of 
140  km.  Though  the  data  is  small,  there  is  no  evidence  for  a  major  low  velocity 
zone  beneath  either  province.  Different  velocity  gradients  between  the  two 
provinces  implies  different  temperature  structure,  which  supports  the  hypothesis 
that  a  deep,  cool  lithospheric  root  exists  beneath  the  Kapvaal-Rhodesian  craton. 
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Report  Summary 


Task  Objectives 

The  general  objective  of  this  research  is  to  understand  the  factors  important  in  shaping  the 
seismic  signature  of  small  events  recorded  at  local  and  regional  distances.  Specific  objectives 
are  1)  the  identification  of  deterministic  aspects  of  the  wavefield  from  small  earthquakes  and 
explosions  to  allow  the  inference  of  source  depth  and  other  source  parameters,  and  2)  to  under¬ 
stand  the  effects  of  scattering  and  lithospheric  heterogeneity  on  the  propagation  of  high  fre¬ 
quency  regional  phases.  The  combined  analysis  of  "deterministic"  and  "stochastic"  wave 
propagation  effects  is  required  to  unravel  the  complexity  of  regional  phases  for  the  purposes  of 
event  discrimination. 

Technical  Problem 

Regional  phases  from  small  events  are  affected  by  complex  interactions  between  source  radi¬ 
ation  and  wave  propagation  effects  due  to  structure  in  the  crust  and  upper  mantle.  Because 
observations  are  confined  to  the  high  frequency  band  (>2  hz),  lithospheric  heterogeneity 
becomes  important  in  shaping  high  frequency  regional  phases.  Typical  wave  lengths  are  much 
shorter  than  the  total  travel  path  and  are  comparable  to  known  geologic  structure.  An  aspect  of 
the  problem  of  regional  wave  propagation  is  examined  here  involving  the  nature  of  coda  asso¬ 
ciated  with  major  arrivals.  Simple  wholespace  scattering  models  are  often  applied  to  the  coda  of 
regional  S  phases  to  deduce  scattering  or  anelastic  attenuation.  Tne  level  and  amplitude  decay  of 
coda  is  a  characteristic  of  the  data  which  seems  to  be  robust  for  particular  regions  and  can  be 
used  to  deduce  source  magnitude,  once  calibrated.  Factors  affecting  the  level  of  scattering  near 
the  source  and  near  the  receiver  are  also  obviously  important  in  yield  estimation  problems  and 
waveform  modeling  for  source  parameters.  It  is  important,  therefore,  to  investigate  the  major 
assumptions  contained  in  these  wholespace  models  and  to  determine  which  are  appropriate. 

Several  problems  have  been  addressed  in  this  research  and  include  analysis  of  P  wave  coda 
from  teleseismic  events,  analysis  of  long-period  Rayleigh  wave  coda  from  explosions  at  NTS, 
and  regional  wave  propagation  of  P  and  S  body  waves.  The  first  two  problems  are  reciprocal  in 
some  respects  and  involve  the  nature  of  scattering  of  waves  in  the  lithosphere  and  how  such  scat¬ 
tering  affects  the  seismogram.  The  analysis  of  body  wave  and  surface  wave  data  using  similar 
formalisms  is  useful  in  characterizing  observed  scattered  coda  waves  behind  the  direct  arrivals 
and  in  investigating  the  degree  of  lithopheric  heterogeneity.  Specific  problems  concerning  tele- 
seismic  P  coda  waves  are  1)  determination  of  the  behavior  of  P  wave  coda,  and  2)  deducing  the 
scattering  mechanisms  involved  in  supporting  high  coda  energy  levels.  Long-period  Rayleigh 
waves  were  studied  to  determine  characteristics  of  coda,  if  any,  and  the  relationship  of  coda  Q 
measurements  to  other  standard  attenuation  measurements.  This  study  is  also  a  prelude  to  deter¬ 
mining  whether  Rayleigh  coda  would  be  a  useful  measure  of  yield.  In  addition  to  studies  of  scat¬ 
tered  waves  several  other  deterministic  studies  of  long  and  short-period  waves  were  performed  to 
determine  vertically  inhomogeneous  earth  structure. 

General  Methodology 

The  general  methodology  of  this  research  involves  theoretical  development  of  scattering 
models  to  apply  to  time  series  data  and  the  simulation  of  wave  propagation  of  in  models  of  plane 
layered  elastic  earth  structure.  In  all  cases  synthetic  seismograms  are  computed  to  compare  with 
the  observed  seismic  data. 

Technical  Results 

The  scattering  of  waves  by  2  and  3D  heterogeneity  in  the  lithosphere  is  a  complex  process 
which  is  poorly  understood.  In  the  first  section  of  this  report,  broadband  teleseismic  P  wave  data 
from  PAS  station  was  studied  to  develop  an  appropriate  scattering  model.  Scattering  under  the 
station  was  seen  to  be  quite  severe  because  of  the  existence  of  a  substantial  tangential  component 
of  motion  comparable  to  vertical  and  radial  motions.  The  level  and  decay  of  coda  was  modeled 


first  using  simple  ID  stochastic  structure  models.  The  results  of  this  study  showed  that  unrea¬ 
sonably  high  ID  velocity  variations,  amounting  to  20%  and  more,  were  needed  to  explain  the 
slow  decay  and  high  level  of  coda.  The  ID  simulations  also  showed  that  the  single  scattering 
model  for  a  simple  plane  wave  in  a  homogeneous  wholespace  was  not  particularly  useful  in 
interpreting  coda  behavior.  In  fact,  the  interpretation  of  coda  decay  with  the  single  scattering 
model  would  imply  that  scattering  attenuation  is  greatest  for  those  structures  which  are  not  heter¬ 
ogeneous.  The  ID  models  showed  that  decay  rate  decreases  with  degree  of  scattering.  Observa¬ 
tional  aspects  of  the  coda  data  were  modeled  using  an  energy  flux  model  for  scattering  in  a  layer 
over  homogeneous  halfspace.  The  energy  flux  model  is  a  useful  parameterization  of  the 
teleseismic  coda  data  which  can  be  used  in  comparative  coda  studies  between  seismic  stations. 

These  unusual  results  prompted  a  study  of  coda  from  regionally  propagating  waves.  The 
teleseismic  scattering  model  suggested  that  apparent  attenuation  deduced  from  coda  decay  is 
strongly  affected  by  the  gross  distribution  of  scatterers  as  well  as  their  density.  Fast  coda  decay 
with  time,  and  consequently  low  Q,  implies  that  waves  quickly  radiate  into  the  mantle  rather  than 
being  attenuated  through  scattering  losses.  Section  2  describes  a  study  of  long-period  Rayleigh 
wave  coda  where  various  mechanisms  of  coda  formation  were  investigated.  There  has  been 
extensive  previous  work  on  attenuation  of  Rayleigh  waves  in  western  North  America  using 
amplitude  decay  of  fundamental  mode  Rayleigh  waves.  The  coda  study  was  attempted  to  first 
determine  if  long-period  Rayleigh  waves  had  coda  (they  do)  and  then  to  investigate  whether  scat¬ 
tering  mechanisms  could  be  separated  from  anelasticity  (they  can’t).  The  single  scattering  model 
and  two  different  energy  flux  models  were  used  to  model  coda.  The  results  show  that  coda 
decay  is  the  most  observationally  robust  measurement  and  that  coda  decay  could  be  due  to  ane- 
lastic  attenuation.  However,  an  energy  flux  model  which  allows  the  scattering  and  subsequent 
radiation  of  surface  wave  to  body  wave  conversions  into  the  mantle  produces  attenuation  effects 
comparable  to  the  scattering  of  teleseismic  P  waves  seen  in  the  PAS  data.  Both  of  these  scatter¬ 
ing  studies  raise  very  interesting  questions  on  the  mechanisms  of  scattering  as  well  as  on  the  role 
of  scatterer  geometry  in  controlling  the  characteristics  of  coda. 

In  conjunction  with  these  studied  on  coda,  an  analysis  of  regional  P  wave  propagation  was 
completed  to  determine  upper  mantle  structure  in  the  South  African  shield  (Section  3).  This  was 
prompted  by  observations  of  anomalous  P^  waveforms  for  crustal  events  believed  to  be  asso¬ 
ciated  with  southern  extension  of  the  east  African  rift  zone  into  Zambia.  Structure  and  wave 
propagation  in  shield  areas  is  a  recurring  problem  in  discrimination  and  detection  since  there  are 
large  shield  areas  of  Asia  and  Europe  which  must  be  monitored.  Past  studies  of  the  P„,  phase 
have  concentrated  in  using  this  phase  to  infer  gross  crustal  structure  and  focal  mechanisms.  In 
the  present  study,  aspects  of  the  observed  PnI  phases  were  used  to  infer  the  character  of  P  wave 
velocity  gradients  in  the  upper  mantle.  The  study  showed  that  substantial  positive  velocity  gradi¬ 
ents  are  required  to  fit  the  large  long-period  P„  waves  observed  in  the  data.  Positive  velocity 
gradients  in  the  uppermost  mantle  give  rise  to  turning  waves  and  high  relative  amplitudes  for 
direct  P  and  S  wave  phases  seen  at  regional  distances.  This  kind  of  structure  is  conducive  for 
detection  of  small  of  events  at  regional  distances.  It  is  also  important  to  recognize  when  such 
structure  exists  between  source  and  receiver  since  magnitudes  and  yields  determined  from  direct 
body  wave  phases  may  be  overestimated  due  to  the  more  efficient  wave  propagation. 

Important  Findings  and  Conclusions 

Analysis  of  coda  in  teleseismic  receiver  functions  suggests  that  there  are  other  mechanisms 
which  control  the  formation  and  decay  of  scattered  coda  waves  which  are  separate  from  intrinsic 
or  scattering  attenuation.  In  particular,  the  simple  geomeoy  of  an  elastic  scattering  layer  over 
halfspace  produces  coda  decay  which  would  be  interpreted  as  an  attenuation  effect  but  is  due  to 
the  simple  redistribution  of  scattered  energy  from  the  layer  to  the  halfspace.  Such  effects  can  be 
studied  first  by  analyzing  coda  from  teleseismic  events  at  a  receiver  or  array  and  then  applying 
the  parameters  of  the  layer  model  to  an  appropriate  energy  flux  model  for  a  source  contained 
within  the  layer. 

Study  of  regional  long-period  surface  propagation  shows  that  several  mechanisms  affect 
coda  decay  and  level.  Although  coda  decay  can  be  explained  and  is  consistent  with  anelastic 


attenuation  measured  previously,  scattering  mechanisms  can  be  assumed  to  create  the  same 
effects.  This  has  important  implications  in  the  interpretation  of  coda  in  all  wave  propagation 
regimes. 

Significant  Hardware  Development 

N/A 

Special  Comments 

N/A 

Implications  for  Further  Research 

Scattering  Q  models  are  based  on  a  number  of  assumptions  concerning  the  distribution  of  the 
reservoir  of  energy  contained  in  the  scattered  field.  Further  research  will  concentrate  on 
combined  application  of  teleseismic  scattering  determiniations  and  regional  phase  scattering. 

This  combined  analysis  may  offer  constaints  on  scattering  physics  not  obtainable  by  analysis  of 
each  data  set  alone.  Results  of  this  research  will  have  important  implications  on  studies  of 
regional  phase  propagation,  discrimination  of  small  events,  and  yield  estimation  problems. 

The  work  performed  on  this  contract  has  generated  several  new  areas  of  research  which  are 
being  pursued.  Figure  1  shows  results  of  a  finite  difference  calculation  designed  to  study  the 
energy  flux  model  developed  here  for  teleseismic  P  waves  and  receiver  scattering.  The  finite 
difference  model  consists  of  a  2D  scattering  layer  over  a  homogeneous  halfspace  (acoustic  case 
only).  The  calculations  are  being  performed  to  first  test  the  assumptions  that  go  into  the  energy 
flux  model  and  then  to  gain  understanding  of  the  nature  of  scattering  in  the  2D  layer  compared  to 
the  ID  structures  used  before.  Simulations  in  solid  elastic  structure  are  also  being  performed  to 
look  at  the  partitioning  of  scattered  waves  into  rotational  and  dilatational  fields  and  how  such 
partitioning  affects  the  interpretations  of  the  parameters  of  the  energy  flux  model.  Such  calcula¬ 
tions  will  also  be  extended  to  scattering  with  a  source  in  the  layer  for  regional  wave  propagation 
studies. 

An  observational  study  of  high  frequency  data  recorded  by  the  NORESS  array  has  also 
been  started  to  investigate  regional  phases  and  the  influence  of  deterministic  (ID  structure)  vs 
stochastic  earth  structure.  Figure  2  displays  several  wave  forms  from  different  events  recorded 
by  NORESS.  The  study  is  taking  advantange  of  the  array  capabilities  by  analyzing  any  individ¬ 
ual  phase  seen  sweeping  across  the  array,  i-k  analysis  anu  beam  forming  is  being  used  to  first 
detect  phases  and  then  to  model  their  characteristics  using  plane  layered  models.  Figure  2  shows 
the  result  of  stacking  individual  time  windows  about  an  observed  phase  at  the  optimal  phase 
velocity  (velocity  show  by  numbers,  time  windows  by  vertical  dotted  lines).  Each  stack  is  then 
patched  together,  for  an  event,  to  display  an  enhanced  seismogram.  Incoherent  scattered  waves 
are  reduced  by  this  method.  The  stacked  data  show  a  number  of  impulsive  phases,  some  with 
quite  high  phase  velocities.  The  high  velocities  can  be  understood  by  considering  the  crusta! 
structure  inferred  for  the  area.  For  example,  the  9.2  km/sec  velocity  is  consistent  with  a  precriti- 
cal  wide  angle  P  reflection  from  the  Moho.  Initial  results  from  this  study  show  that  a  number  of 
unusual  phases  occur  in  the  data  associated  with  multiples  and  conversions  within  the  crust. 
These  phases  can  be  used  to  locate  the  event  as  well  as  determine  crustal  structure.  Ultimately, 
we  want  to  be  able  to  understand  the  short-period  regional  seismogram  to  find  near-source  depth 
phases  for  discrimination  purposes.  This  particular  study  is  also  a  prelude  to  a  larger  study  of 
scattering  of  teleseismic  and  regional  P  waves  at  NORESS,  following  up  on  several  issues  asso¬ 
ciated  with  the  energy  flux  model  developed  in  this  contract. 


Figure  1  -  Layer  over  halfspace  model  (try)  of  a  scattering  layer  with  a  random  velocity  described  by  a  10%  standard  deviauon  Averge  P  a  e 
velocities  for  the  laver  .nd  Mfspec-  ,-T,nd  *  lrm/«-c  respectively.  The  velocity  scale  shows  deviations  from  the  average.  The  lower  dime 
traces  are  synthetic'seismograms  for  receivers  at  the  surface  near  the  center  of  the  model.  A  vertically  incident  P  wave  was  assumes,  .  ole  the 
long  duration  of  coda  and  its  near-constant  level. 
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Stacked  Vertical  Components 


Figure  2  ~  Enhanced  vertical  seismograms,  normalized  to  maximum  amplitude 
plotted  as  a  function  of  distance.  The  stacking  velocity  within  each 
slant-stack  time  window  is  indicated. 
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Scattering  of  Teleseismic  Body  Waves  Under  Pasadena,  California 

Charles  A.  Langston 
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Teleseismic  receiver  funcnons  For  *"-uciure  under  PasaJena.  California  iPAS)  are  derived  from  azi- 
muihally  distributed  teleseismic  P  waves  recorded  on  BenioiF  1-90  instrumentation.  The  broadband 
three-component  Benioff  1-90  system  is  peaited  at  a  1-s  period  and  allows  resolution  of  major  crustal 
interfaces  from  large  Ps  conversions  seen  in  the  receiver  function  data  The  observed  body  wave  data  are 
quite  complex,  showing  exceptionally  large  Ps  conversions  and  scattered  waves  on  horizontal  compo¬ 
nents  Radial  and  tangential  motions  are  of  equal  magnitude  and  show  major  olf-azimuth  converted  Ps 
wases.  suggesting  large-scale  crustal  heterogeneity  beneath  the  station  Stochastic  simulations  of  one- 
dimensional  plane  layered  structure  show  that  geologically  unreasonable  one-dimensional  models  are 
required  to  fit  the  data.  The  observed  coda  decay  yields  a  scattering  Q  estimate  of  239  at  a  2-s  penod 
using  an  energy  flux  model  for  a  propagating  plane  wave  interacting  with  a  scattering  layer  over  a 
homogeneous  half-space  Observed  and  synthetic  coda  decay  follows  the  theoretical  exponential  decay 
predicted  by  the  model  and  is  due  entirely  to  diffusion  of  coda  energy  out  of  the  layer  into  the  half-space 
PAS  coda  is  compared  to  coda  from  deep  teleseisms  recorded  at  State  College.  Pennsylvania,  and  it  is 
seen  that  scattering  is  more  severe  at  PAS,  as  reflected  in  higher  coda  levels  and  slower  decay  rate 
Consideration  of  energy  pariibomng  and  coda  amplitude  suggests  that  much  of  the  coda  consists  of 
scattered  surface  waves  Analysis  of  a  major  Ps  conversion  arriving  3  s  after  direct  P  indicates  that  a 
major  crustal  discontinuity  at  about  20  km  depth  dips  at  moderate  angles  to  the  north  under  the  San 
Gabriel  Mountains.  This  interface  probably  represents  the  crustal  tectonic  boundary  between  the  Trans¬ 
verse  Ranges  and  the  Los  Angeles  Basin 


Introduction 

The  analysis  of  teleseismic  receiver  functions  represents  an 
inexpensive  and  convenient  way  of  imaging  major  crustal  and 
upper  mantle  discontinuities  under  isolated  receivers.  The 
transmissivity  of  structure  under  a  three-component  seismo¬ 
meter  is. inferred  from  the  timing  and  amplitude  of  Ps  conver¬ 
sions  seen  on  horizontal  ground  motions  and  is  modeled  to 
determine  the  location  and  velocity  contrasts  of  the  causitive 
interfaces  [ Burdtek  and  Langston,  1977;  Langston.  1979; 
Owens.  1984],  The  technique  has  been  particularly  useful  in 
large-scale  structure  studies  using  long-period  body  waves 
[eg..  Burdick  and  Langston.  1977;  Langston  and  Isaacs.  1981; 
Hebert  and  Langston.  1985]  to  determine  average  crustal 
thickness  and  is  increasingly  being  applied  to  broadband, 
high-frequency  data  to  obtain  more  resolution  on  structure 
[Owens.  1984;  Owens  c(  at..  1984;  1987] 

One  of  the  inevitable  trade-offs  in  using  higher-frequency 
data  is  increased  sensitivity  to  lateral  heterogeneity  in  crustal 
structure  In  one  sense  this  is  desirable  since  a  goal  of  such 
studies  is  to  determine  as  much  information  about  structure 
under  the  receiver  as  possible.  However,  it  is  also  obvious  that 
the  wave  field  is  severely  spatially  aliased  through  observa¬ 
tions  made  at  only  a  single  surface  point  Imaging  procedures 
implicitly  rely  on  modeling  assumptions  such  as  plane  layering 
or,  at  most,  simple  curved  interfaces.  It  i;  often  observed  that 
much  of  the  wave  field  is  inaccessable  to  standard  explanation 
..sing  simple  modeling  techniques  [Langston.  1979;  Owens  er 
u/.,  1987]  For  example,  receiver  function  data  often  display 
anomalous  wave  behavior  such  as  P  wave  particle  motions 
which  have  significant  tangential  amplitudes. 

A  purpose  of  this  paper  is  to  examine  strategies  of  treatment 
of  broadband  receiver  function  data  which  take  into  account 
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both  stochastic  and  deterministic  scattering  effects  due  to  het¬ 
erogeneous  structure  Structure  under  the  station  Pasadena. 
California  I  PAS)  will  be  the  focus  of  this  effort.  This  station 
lies  in  a  geologically  complex  area  of  southern  California  and 
has  had  broadband  instrumentation  over  a  long  period  of 
time  (Figure  If  The  study  was  initially  motivated  by  observa¬ 
tions  of  complex  particle  motions  in  the  teleseismic  P  wave 
data  [Langston.  1977],  The  data  were  previously  examined  by 
Lee  [1983],  who  attempted  to  model  the  receiver  functions 
using  dynamic  ray  tracing  with  models  consisting  of  homoge¬ 
neous  layers  separated  by  curved  three-dimensional  interfaces 
Lee  s  study  was  partially  successful  in  explaining  qualitative 
aspects  of  the  data  for  early  arrivals,  but  he  found  that  ray 
theory  was  inadequate  to  explain  the  high  amplitude  of  in¬ 
ferred  covcrted  waves  and  the  duration  of  signal. 

In  this  paper  the  receiver  function  data  will  be  examined 
from  two  points  of  view  The  first  is  from  the  standard  method 
of  treating  the  data  to  infer  major  velocity  discontinuities 
under  the  station  using  observations  of  isolated  body  wave 
phases  and  simple  velocity  models.  Ps  conversions  from  telc- 
seismic  P  waves  are  found  to  be  unusally  large  and  are  used  to 
suggest  the  existence  of  a  large  velocity  contrast  interface  in 
the  lower  crust  which  dips  to  the  north  under  the  San  Gabriel 
Mountains 

The  other  point  of  view  is  to  treat  the  data  as  resulting  from 
an  unknown  scattering  process  and  to  attempt  to  infer  the 
severity  of  wave  scattering  under  the  station  Using  simple 
measures  of  the  P  coda  amplitude  decay  along  with  one- 
dimensional  stochastic  structure  simulations,  the  question  is 
asked  Are  the  data  consistent  with  scattering  due  to  reason¬ 
able  plane  layered  structure  ’  The  negative  answer  for  PAS 
suggests  that  such  an  analysis  can  be  used  to  justify  or  not 
justify  a  research  effort  in  modeling  data  with  simple  plane 
layered  structure  models  The  severity  of  observed  scattering 
under  PAS  also  points  out  deficiencies  in  some  simple  scatter¬ 
ing  models  and  the  need  to  develop  appropriate  models  for 
two-  and  three-dimensional  stochastic  structures. 
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Fig  I  Sketch  map  of  southern  California  showing  the  lc-  itton  of 
PAS  station  and  major  faults  of  the  area 


PAS  Station  and  Regional  Structure 

PAS  station  has  been  in  operation  since  the  mid- 1 930s  and 
has  had  a  full  complement  of  experimental  long-  and  short- 
period  instrumentation.  Of  interest  to  this  study  are  ihe  Be- 
moff  1-90  (seismometer  period  is  1  s.  galvonometer  period  is 
90  s)  and.  to  a  lesser  extent,  the  Press-Ewing  30-90  systems. 
The  Benioff  1-90  system  is  peaked  at  1-s  period  but  records 
across  a  wide  frequency  band  comparable  to  the  intermediate- 
period  Digital  World-Wide  Standard  Seismograph  Network 
(DWWSSN)  passband  (see  the  appendix).  The  system  has  rou¬ 
tinely  recorded  teleseisms  throughout  its  existence  Although 
the  data  are  recorded  in  analog  photographic  format,  the 
broad  passband  of  the  instrument  allows  for  significant  time 
resolution  of  crustal  Ps  conversions.  This  potential,  in  con¬ 
junction  with  the  complexity  of  the  receiver  signal  [ Langston , 
1 977 J  as  well  as  the  tectonic  problems  associated  with  the 
Transverse  Ranges  Province  of  southern  California,  motivates 
the  present  study  of  crustal  structure  under  the  receiver. 

PAS  station  lies  near  the  southern  boundary  of  ihe  Trans¬ 
verse  Ranges  Province  and  the  Los  Angeles  Basin- Peninsular 
Range  Province  (Figure  I).  Crustal  thickness  from  long-range 
refraction  [ Hadley  and  Kanamori,  1977]  and  time  term  analy¬ 
sis  [ Hearn  and  Clayton,  1986h]  suggest  that  the  crust  under 
PAS  is  about  31  km  thick.  Geologic  structure  of  the  upper 
crust  is  known  to  be  quite  complex  with  major  active  faults 
separating  regions  of  diverse  rock  type.  For  example,  crys¬ 
talline  rocks  of  the  San  Gabriel  Mountains  a  few  kilometers 
north  of  PAS  abut  valley  sediments  that  attain  depths  of  up  to 
10  km  in  the  surrounding  basins  [Writs  et  al..  1985],  The 
geology  of  the  Transverse  Ranges  and  the  San  Gabriel  Moun¬ 
tains  suggest  that  these  ranges  are  in  part  allochthonous, 
being  thrust  over  younger  rocks  of  the  Peninsular  Range 
Province. 

Hadley  and  Kanamon  [1977]  review  a  number  of  long- 
range  refraction  and  travel  time  studies  for  the  area  and  show 
that  the  Transverse  Ranges  are  a  locus  of  both  upper  mantle 
and  crustal  velocity  anomalies.  Using  the  Southern  California 
Seismic  Network,  they  showed  that  P  delays  from  a  PKIKP 
ph  *  outlined  an  east-west  trending  /one  of  high  velocities  in 


Ihe  upper  mantle  This  /one  was  studied  hy  Humphre  i  v  et  al 
[1984]  using  a  lomographc  imaging  technique  for  teleseismie 
P  waves  recorded  h\  the  network  They  determined  that  the 
vertical  zone  ol  high  mantle  velocities  wa-  contained  within 
the  Transverse  Ranges  Province  and  attained  depths  of  at 
least  150  km  This  zone  was  interpreted  as  a  region  of  mantle 
downwellmg  analogous  to  a  subduction  zone  but  driven  by 
relative  movements  of  microplates  making  up  the  San  An¬ 
dreas  fault  system.  In  this  model,  upper  crustal  rmcroplate 
movements  may  be  decoupled  from  lower  crust  and  upper 
mantle  structures  with  the  result  that  the  surface  expression  of 
the  San  Andreas  fault  is  offset  by  a  midcrustal  horizontal 
shear  zone  from  the  fault  at  depth  [ Webb  and  Kanamon. 
1985] 

The  Transverse  Ranges  are  3lso  a  locus  of  change  in  crustal 
structure  between  the  western  Peninsular  Ranges- Los  Angeles 
Basin  and  the  Mojave  Block  to  the  east.  Hadley  and  Kanamon 
[1977]  suggest  that  a  high-velocity  lower  crustal  layer  com¬ 
prises  about  half  of  the  crust  in  the  Peninsular  Ranges  but 
tapers  to  onlv  a  few  kilometers  in  the  Mojave  Block  Tomo¬ 
graphic  study  of  Pa  and  Pn  waves  m  southern  California 
[Hearn  and  Clayton ,  1986u.  h ]  support  this  suggestion  by 
showing  slower  average  velocities  in  the  Mojave  Block  relative 
to  crust  to  the  west. 

A  major  goal  of  the  present  work  was  initially  to  provide 
constraints  on  crustal  thickness  un-er  this  important  transi¬ 
tion  between  tectonic  provinces  Site  specific  information  pro¬ 
vided  by  the  receiver  function  technique  will  complement 
these  earlier  crustal  and  upper  mantle  structure  studies  but 
will  also  show  the  presence  of  a  major  midcrustal  to  lower 
crustal  interface. 

Data  and  Sourcf  Function  Equalization 

Wave  form  data  from  21  teleseismie  earthquakes  were  ob¬ 
tained  from  the  seismogram  archives  of  California  Institute  of 
Technology  (Table  1).  The  seismograms  were  photographed 
ar.i  enlarged  for  hand  digitization.  Wave  forms  were  digitized 
at  an  irregular  sampling  interval  and  interpolated  to  an  equal 
sampling  interval  of  0.1  s  Processing  included  vector  rotation 
of  the  horizontal  components  into  the  theoretical  back  azi¬ 
muth  of  the  P  wave  arrival  to  obtain  radial  l positive  away 
from  the  sourcel  and  tangential  (positive  clockwise  around  the 
source  looking  at  the  receiveri  ground  motions 

A  source  function  equalization  procedure  was  ihen  per¬ 
formed  to  remove  the  instrument  respo-  .e  and  unknown  ef¬ 
fective  source  function  from  the  m  .al  and  tangential  wave 
form  data  [Langston.  1979]  this  procedure  the  vertical 
component  of  motion  is  .ssumed  to  be  free  of  any  effect  of 
near-receiver  structure  (reverberations  or  conversions)  but 
contains  the  common  instrument  response  and  wave  propaga¬ 
tion  effects  from  the  mentle  and  near-source  region  The  data 
are  time-windowed  and  Fourier-transformed  The  vertical 
spectrum  is  divided  into  the  horizontal  spectra  and  then 
multiplied  by  a  Gaussian  function  to  remove  high-frequenev 
noise.  The  spectral  division  is  also  accompanied  by  prewhiten- 
ing  the  vertical  component  spectra  using  a  "water  level1'  pa¬ 
rameter  to  remove  spurious  spectral  holes  The  water  level 
used  for  data  considered  here  was  0.1".,  of  the  maximum  of 
the  vertical  component  amplitude  vpectrum  The  Gaussian 
filter  used  is  equivalent  to  a  Gaussian  pulse  in  the  time 
domain  with  a  half  width  of  I  s  (i  e  .  u  =  I  f  in  e‘J"',:l. 

F  igure  2  shows  examples  of  wave  form  data  for  events  in 
the  three  major  back  azimuth  i-inges  128  .  235  .  and  315 
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TABLE  I.  Fvenl  Parameters  for  Pasadena  Data 


Origin  Dale 

Time. 

LIT 

Latitude. 

deg 

Longitude. 

deg 

w* 

De  pi  h . 
km 

Distance, 

deg 

Back  Azimuth, 
deg 

Stack 

Group 

March  17.  1966 

1550  33 

21  IS 

179  2  W 

6.2 

639 

79  V 

236.1 

235 

Nov.  20.  1971 

0728:01 

13.4S 

179  9W 

5M 

81.9 

234  8 

235 

Feb.  1.  1973 

0514  20 

22  7S 

62  2W 

6  1 

22V 

77  6 

128  4 

128* 

Dec  28.  1973 

0531  06 

23. 9S 

180  OK 

6.3 

5-19 

82  3 

234.5 

235- 

March  23.  1974 

1428:35 

23  9S 

179  8E 

6.1 

535 

82.5 

234  6 

235* 

Oct.  20.  1974 

0412:29 

17. 9S 

1 78.635' 

6  0 

602 

77.4 

238.1 

235 

Nov.  29.  1974 

2205:22 

30. 7N 

I38.3F. 

6  1 

419 

83  3 

302.5 

ID* 

Feb.  22.  1975 

2204:37 

24. 9S 

179  IW 

6.2 

375 

82.4 

233.7 

235 

June  29.  1975 

1037:4) 

38. 8N 

6  2 

560 

83.9 

313.2 

315 

Jan.  23.  1976 

0545:30 

7.5S 

1I9  9E 

fv4 

614 

120.0 

282.4 

235 

Dec.  12.  1976 

0108:50 

28. ON 

I39.6F 

5.9 

491 

83  9 

299.7 

315 

Sept.  23.  1978 

1644:26 

11.  OS 

167  2E 

6  5 

200 

85.4 

249  9 

235 

April  24,  1979 

0156.14 

20  8S 

178.7W 

ft.O 

450 

79.4 

236.0 

235 

May  13.  1979 

0638:15 

I8.9N 

145. 3E 

5  9 

250 

84.8 

289.2 

315 

Mav  21.  1979 

2232:58 

15  2S 

70. 1W 

M) 

208 

67,1 

128.8 

128 

June  27.  19'9 

0958:03 

7. IN 

82. 0W 

5.8 

150 

42.9 

120.6 

128 

Aug  16.  1979 

2142:44 

41  8N 

(30.7E 

6  I 

588 

81.5 

315.2 

315 
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Also  shown  are  PKIKP  wave  forms  of  the  January  23,  1976. 
event  which  was  used  in  Hadley  and  Kanamori's  [1977]  study. 
This  phase  will  be  used  as  a  contraint  on  structure  under  the 
station  Note  that,  in  all  cases,  the  magnitude  of  tangential 
motions  excited  by  the  P  wave  is  comparable  to  radial  mo¬ 
tions  In  a  radially  stratified  earth.  P  and  resulting  P-SF  con¬ 
versions  would  be  restricted  to  the  sagittal  plane  containing 
the  ray  Note  also  that  tangential  and  radial  motions  differ 
quite  substantially  in  wave  form,  suggesting  that  simple  instru¬ 
ment  miscalibration  or  magnification  differences  cannot  give 
tse  to  these  anomalous  particle  motions  (see  the  appendix). 

The  PKIKP  phase  fo  the  January  23,  1976,  event  also 
shows  very  anomalous  particle  motions.  This  phase  is  incident 
below  the  crust  at  an  incident  angle  of  about  4°.  Horizontal 
motions,  however,  are  not  small  as  expected.  They  are  about 
half  the  size  of  the  vertical  P  wave,  and  both  components  are 
grossly  different  Furthermore,  this  event  has  a  back  azimuth 
in  which  the  observed  E-W  component  is  almost  perfectly 
radial,  and  the  observed  N-S  component,  tangential.  The  data 
of  Figure  2  strongly  suggest  that  heterogeneous  three- 
dimensional  structure  is  causing  large  scattering  effects  in  the 
receiver  function  data. 

Because  of  the  location  of  teleseismic  source  zones,  the 
equalized  radial  and  tangential  component  data  were  grouped 
into  three  back  azimuth  groups  and  stacked  (see  Table  1  for 
groupings).  The  stacking  was  done  by  shifting  all  traces  to  a 
common  time  based  on  the  P  first  arrival  and  then  adding 
them  Wave  forms  for  one  standard  deviation  about  the  mean 
at  each  sample  point  were  also  computed.  The  resulting  wave 
form  stacks  are  shown  in  Figure  3  with  their  standard  devi¬ 
ations  A  comparison  of  stacks  from  the  three  back  azimuth 
groups  are  shown  in  Figure  4.  Displays  of  this  type  yield 
information  on  the  coherency  of  arrivals  within  the  wave  form 
and  of  the  level  of  processing  noise  [Owens,  1984] 

Only  the  first  15  s  of  the  wave  forms  are  shown  in  F  igure  3. 
The  data  of  Figure  2  show  major  arrivals  in  the  horizontal 
wave  forms  for  at  least  60  s  This  poses  a  dilemma  for  velocity 
models  that  can  be  considered,  since  it  is  difficult  to  gel  such 
arrivals  from  plausible  plane  layered  models  This  problem 


will  be  addressed  in  a  later  section;  here  I  concentrate  on 
major  initial  arrivals 

Figure  3 ci  shows  the  wave  form  stacks  for  events  from  a 
back  azimuth  of  128  The  bounding  envelope  for  one  stan¬ 
dard  deviation  is  large  for  tangential  motions  and  for  arrivals 
after  direct  P  on  the  radial  component  A  phase  marked  "Ps" 
on  the  radial  stack  is  also  observed  on  the  radial  data  of  the 
other  back  azimuth  groups. 

The  other  two  back  azimuth  groups  (Figures  3 h  and  3c) 
show  remarkably  large  arrivals.  The  Ps  conversion  is  roughly 
half  the  size  of  direct  P  on  the  radial  components  and  is  also 
large  on  the  tangential  components.  Note  that  the  tangential 
components  show  that  the  direct  P  wave  amplitude  is  variable 
within  the  noise  of  measurement.  There  are  indicatio  s  of 
later,  coherent  arrivals  in  the  wave  forms,  but  these  are  f  ib- 
lematical 

Figure  4  shows  the  stacks  displayed  together  with  P  ana  Is 
phases  annotated  My  working  hypothesis  is  that  this  secor  ' 
ary  phase  is  a  direct  P-to-SF  conversion  beneath  the  statio 
Its  large  size,  relative  to  direct  P.  on  all  except  the  128 
tangential  stack  is  remarkable  and  can  be  directly  seen  in  the 
data  of  Figure  2  (also  see  the  appendix).  Figure  5  shows  parti¬ 
cle  motion  plots  of  the  radial  and  tangential  stacks  for  315 
and  235°  back  azimuths.  Note  that  the  direct  P  wave  con¬ 
forms  to  nearly  radial  motions  as  expected  for  ideal  P  particle 
motion  but  that  the  Ps  conversion  has  been  rotated  45  or 
more  out  of  the  sagittal  plane.  It  is  very  difficult  to  produce 
such  Ps  arrivals  from  simple  dipping  interfaces  that  dip  only  a 
few  degrees  or  have  velocity  contrasts  typical  of  continental 
crustal  models  [Langston,  1977;  Lee.  1983].  The  implication  is 
that  there  is  a  major  discontinuity  under  PAS  which  has  high 
dip  and/or  high  S  wave  velocity  contrast. 

These  arrivals  are  also  evident  in  long-period  data.  Figure  6 
compares  observed  wave  forms  for  the  November  29,  1974. 
event  recorded  on  (he  BemofT  1-90  and  Press-Twang  30-90 
systems  at  PAS  The  data  have  been  shifted  10  a  common  lime 
base  The  vertical  components,  although  showing  some  scat 
tered  waves  in  ihc  coda,  are  pulselike  and  show  one  major  P 
arrival  The  N-S  1-90  component  is  dominated  by  the  Pi  con- 
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Fig  2  Selected  data  (rum  deep  earthquakes  (see  Table  ll  recorded  on  BemolT  1-90  insirumemalior  P6S  7.  N  S. 
and  E-W  wave  forms  denoie  observed  vertical,  norih-soulh.  and  east-west  components  R  and  T  «?n  ,urms  are  the  result 
of  rotating  the  observed  horizontal  components  into  the  theoretical  back  azimuth  of  the  rav  I  he  distance  and  back 
azimuth  angle  arc  given,  in  order,  in  the  parenlheses  to  the  right  of  each  event's  dale  Note  ;ne  lime  scale  difference  for  the 
January  23.  1976.  event. 


version,  and  direct  P  is  a  minor  initial  arrival.  The  30-90  N-S 
record  shows  that  the  Ps  conversion  (arrowl  also  dominates 
the  long-period  wave  form.  The  Ps  conversion  is  present  on 
the  E-W  components,  hut  direct  P.  the  first  pulse,  is  larger. 
The  Ps  conversion  is  evident,  having  the  effect  of  broadening 
the  initial  pulse  by  a  factor  of  2  compared  to  the  vertical 
long-period  P  wave,  producing  a  shoulder  on  this  pulse  This 
comparison  of  data  recorded  on  two  different  seismometer 
systems  demonstrates  that  the  crustal  structure  responsible  for 
these  scattering  effects  is  large  and  that  the  effect  is  not  an 
artifact  of  instrument  miscalibration. 

Langston  [1977,  1979]  suggests  that  the  amplitude  and  po¬ 
larity  behavior  with  azimuth  of  tangential  Ps  can  constrain 
the  direction  and  magnitude  of  interface  dip.  Unfortunately. 


this  phase  is  only  well  developed  in  the  235  and  315  azimuth 
groups,  although  it  may  be  significant  that  it  is  poorly  devel¬ 
oped  in  the  128  azimuth  group.  The  PKIKP  phase  from  the 
January  23  1976  event,  however,  offers  some  independent 
information  in  this  regard 

Since  this  phase  propagates  nearly  vertically,  any  Ps  con¬ 
version  from  a  dipping  interface  will  be  contained  in  the  plane 
of  the  ray  and  dip  direction  Thus  a  simple  plot  of  particle 
motion  of  the  Ps  conversion  will  be  a  direct  measurement  of 
the  dip  direction  (FTgure  71  The  Ps  conversion  is  the  largest 
arrival  on  the  N-S  component  and  is  easily  seen  in  the  particle 
motion  plot.  Note  that  il  is  polarized  almost  perfectly  north¬ 
ward  Because  Ps  is  the  same  polarity  as  direct  radial  P 
(Figure  4),  it  is  due  to  conversion  at  the  boundary  between 
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Fig  3  Stacks  of  source-equalized  radial  Hop)  and  tangential  ibottom)  data  for  the  three  azimuthal  groups  considered 
in  the  test.  The  average  and  ±  1  standard  deviation  wave  forms  are  shown  in  each  panel  P  and  Ps  arrivals  are  annotated, 
lal  128  group  (b)  235  group  (cl  315  group. 


deeper  high  velocity  material  under  shallower  lower  velocity 
material.  Taking  the  negative  polarity  of  the  PKIKP  phase 
into  account  yields  a  northward  dip  direction  for  the  postu¬ 
lated  interface.  This  is  consistent  with  the  tangential  Ps  polari¬ 
ties  displayed  by  the  235’  and  315°  stacks. 

Deterministic  Modeling  of  the  Ps  Conversion 

Recent  efforts  in  modeling  receiver  function  data  have  con¬ 
centrated  on  formal  inversion  of  the  data  to  obtain  plane¬ 
layered  crustal  and  upper  mantle  models  [e.g.,  Owens  et  ai, 
1987],  Characteristics  of  the  Pasadena  data  set  preclude  this 
approach.  The  duration  and  amplitude  of  the  horizontal  com¬ 
ponent  coda  and  large  tangential  amplitudes  all  argue  against 
finding  reasonable  plane  layered  models.  Figure  8  shows  a 
comparison  of  observed  and  synthetic  radial  component  wave 
forms.  The  radial  stack  for  235°  azimuth  is  shown  below  a 
synthetic  radial  seismogram  computed  for  a  crustal  model 
proposed  for  the  area  by  Hadley  and  Kanamori  [1977],  The 
Thompson-Haskell  method  was  used  to  construct  the  synthet¬ 
ic  [Haskell,  1962]  The  crustal  model  is  shown  in  Figure  9 
with  parameters  tabulated  in  Table  2.  The  Moho  occurs  at  31 
km  depth  and  produces  a  moderately  large  Ps  conversion 
which  arrives  4  s  after  direct  P  (shown  by  arrow  in  Figure  8). 
The  observed  Ps  conversion  is  larger  and  arrives  at  least  1  s 
earlier,  suggesting  several  other  models. 

First,  if  tne  Ps  conversion  is  from  the  Moho.  then  the  crust 
must  be  at  most  27  km  thick  if  Hadley  and  Kanamori's  veloci¬ 
ties  are  assumed  However.  Hearn  and  Clayton's  [I986u]  study 
using  Pi/  waves  suggested  average  crustal  P  wave  velocities  in 


the  area  of  about  6.3  km/s  and  a  crustal  thickness  of  about  31 
km.  It  is  possible  that  the  receiver  function  data  are  sampling 
a  local  anomaly  unresolved  by  Hearn  and  Clayton's  data. 
Alternatively,  if  the  crust  is  31  km  thick,  then  the  average  S 
wave  velocity  in  the  crust  must  be  at  least  4  km  s.  The  corre¬ 
sponding  P  wave  velocity  would  be  high  at  6.9  km/s.  assuming 
Poisson  ratios  near  0.25  appropriate  for  crustal  rocks.  In 
either  the  case  of  a  slow  thin  crust  or  a  thicker  high-velocity 
crust  there  should  be  an  anomalous  mass  excess  in  the  crust 
and  upper  mantle  column  which  would  show  up  in  the  gravity 
field.  No  such  anomaly  is  observed  [ Hadley  and  Kanamori. 
1977], 

One  simple  solution  to  this  problem  is  to  accept  the  thick¬ 
ness  and  crustal  velocities  from  previous  studies  and  to  infer 
that  a  mtdcrustal  interface  causes  the  large  Ps  conversion.  A 
plane  layered  model  which  shows  this  hypothesis  is  plotted  in 
Figure  9  land  Table  2)  with  the  corresponding  radial  synthetic 
seismogram  in  Figure  8  Ps  arrivals  from  the  Moho  are  mini¬ 
mized  by  making  the  structure  approximate  a  smooth  gradi¬ 
ent  in  that  region.  The  large  Ps  amplitudes  observed  require  a 
high  S  wave  velocity  contrast  This,  in  turn,  implies  a  velocity 
inversion  in  the  midcrust  The  synthetic  seismogram  shown  in 
Figure  8  for  this  kind  of  low-velocity  zone  (LVZ1  structure 
approximates  the  arrival  time  and  the  double-peak  character 
of  the  actual  data. 

Obviously,  the  plane  layered  model  does  not  explain  the 
anomalous  particle  motion  of  the  Ps  conversion  A  series  of 
ray  theory  calculations  were  performed  to  test  the  dipping 
interlace  model  [ l.anaston,  1977],  The  tangential  Ps  data  re- 
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quire  that  the  interface  dip  northward  under  the  San  Gabriel 
Mountains.  Experience  with  such  calculations  indicates  that 
there  can  be  considerable  trade-off  between  interface  dip  and 
velocity  contrast  [eg.  Langston.  1979],  Waves  which  ap¬ 
proach  a  dipping  interface  from  the  downdip  direction  will 
have  an  effective  angle  of  incidence  which  is  larger  than  waves 
approaching  a  horizontal  interface  This  will  produce  a  larger 
conversion.  Waves  approaching  from  the  updip  direction  will 
tend  to  have  lower  angles  of  incidence  with  less  conversion. 

A  number  of  dipping  interface  models  based  on  the  Hadley 
and  Kanamori  crustal  model  of  Table  2  were  examined.  The 
top  of  the  6.8  km  s  layer  was  allowed  to  dip  up  to  40  Two 
rays  were  traced  through  the  model  These  were  direct  P  and 
the  Ps  conversion  from  the  dipping  interface.  It  was  quickly 
seen  that,  although  it  may  be  possible  to  produce  large  Ps 
conversions  for  rays  which  approach  the  structure  from  the 
downdip  direction,  models  with  dips  greater  than  10  consis¬ 
tently  produced  low-amplitude  Ps  conversions  for  rays  trav¬ 
eling  from  the  updip  direction.  Indeed,  for  P  velocity  contrasts 
of  6. 2-6.8  km/S  and  5.0-6. 8  km  s  (assuming  a  Poisson  solid), 
dips  of  30  resulted  in  Ps  conversions  which  had  opposite 
polarities  relative  to  direct  radial  P.  inconsistent  with  the  data 
(Figure  4)  Thus  interface  dip  should  be  of  the  order  of  10  or 
less.  The  5  wave  velocity  contrast  should  also  be  greater  than 
I  km  s.  The  calculated  Ps  P  ratio  for  updip  ray  incidence  is 
0.19  for  the  5  0-6.8  km  s  interface  with  10  dip  The  observed 
Ps  P  ratio  for  the  315  stack  is  0.57  (Figure  4),  which  represent 
waves  coming  updip  but  at  an  angle  from  the  northward  dtp 
direction.  Calculated  tangential  amplitudes  for  the  Ps  conver¬ 
sion  are  comparable  to  the  radial  amplitudes  and  agree  with 
the  45  polarization  anomaly  seen  in  the  data 


Stochastic  Sixcctcre  Modk  inc, 

Up  to  this  point  the  data  have  been  ireated  from  a  deter¬ 
ministic  point  of  view  The  data  show  thjt  the  first  Ps  conver¬ 
sion  is  indeed  a  major  wave  propagation  ellcu.  and  n  is  rea¬ 
sonable  to  assume  that  other  early  arrivals  will  he  due  to 
direct  interactions  with  discontinuities  in  the  structure  How¬ 
ever,  the  Ps  conversion  is  simply  the  first  of  many  large  arriv¬ 
als  in  the  horizontal  P  wave  coda  Are  these  later  arrivals 
fundamentally  different  from  early  arrivals  '  What  do  these 
large  arrivals  imply  about  the  heterogeneity  of  the  structure 
under  the  station'1  In  particular,  since  there  is  recent  emphasis 
on  modeling  such  data  with  formal  inversion  methods,  can 
plane  layered  structure  mimic  the  coda  seen  in  the  data  and.  if 
so.  what  are  the  implications? 

First,  I  make  the  assumption  that  all  arrivals  seen  after 
direct  P  in  the  data  represent  waves  scattered  in  structure  near 
the  receiver.  This  assumption  is  probably  poor  for  shallow- 
earthquake  sources  simply  because  of  known  propagation  ef¬ 
fects  like  near-source  surface  reflections.  Deep  events,  how¬ 
ever.  should  be  less  affected  by  near-source  scattering. 

For  deep  teleseisms,  tangential  amplitudes  in  the  coda  are 
as  large  or  larger  than  the  vertical  coat  Therefore  if  these 
coda  waves  are  due  to  scattering  |P  to  P  or  S  to  P t  in  struc¬ 
ture  near  the  source  they  must  finally  convert  to  P  waves  '.o 
arrive  soon  after  the  direct  P  at  the  receiver.  If  they  are  P 
waves,  they  must  have  azimuth  anomalies  greater  than  45 
since  they  are  so  large  on  (he  tangential  component  This 
contradicts  the  near-source  scattering  hypothesis  since  the 
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Fig  4  Comparison  of  the  stacked  radial  and  tangential  equalized 
wave  forms  for  the  three  hack  azimuth  groups  Note  the  azimuthal 
dependence  of  the  Ps  conversion  on  the  radial  components  The  wave 
forms  have  been  shifted  in  baseline  for  viewing  purposes 
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Fig  5.  Panicle  motion  plots  for  the  315  and  235  equalized  wave  forms  The  wave  form  data  are  displayed  above  the 
radial  iR  or  R ADlUangential  (T  or  TANG)  particle  motions.  Arrows  are  shown  on  the  particle  motion  plots  every  0.5  s 
Wave  form  data  included  within  the  brackets  are  plotted  below.  Maximum  amplitude  for  each  wave  form  is  shown  (o 
nght  of  the  wave  form.  Note  that  the  P  waves  in  both  cases  are  polarized  in  the  expected  ray  direction  but  the  inferred  Ps 
conversions  have  a  polarization  anomaly  of  45  . 


scattered  waves  themselves  must  originate  at  teleseismic  dis-  normally  distributed  velocity  values  r|r)  with  zero  mean  and 

tances  from  the  source.  The  conclusion  is  that  large  horizontal  unit  variance  was  generated  from  a  pseudo  random  number 

coda  amplitudes  relative  to  horizontal  P  must  be  due  to  scat-  generator.  A  0  25  km  sampling  interval  was  assumed  for  a 

tertng  near  the  receiver  (see  also  Cessaro  and  Butler  [1987]).  layer  of  thickness  30  or  60  km.  The  velocity  function  was 
Synthetic  seismograms  were  computed  for  a  series  of  receiv-  Fourier-transformed  to  the  wave  number  domain  to  obtain 

er  models  with  plane  layered  stochastic  structure.  The  pro-  iik).  An  exponential  correlation  function  N(r)  =  e-*"  was  as- 

cedure  used  by  Frankel  and  Clayton  [1986]  was  adopted  to  sumed  for  the  medium  where  a  is  the  correlation  length.  The 

generate  a  random  velocity-depth  function.  A  random  series  of  wave  number  spectrum  of  ;V(z),  5?(k)  =  2a/(l  +  k2a2),  was  then 
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Fig  6  Comparison  of  data  from  the  November  29.  1974.  event  recorded  on  30-90  (long  period)  Oop)  and  1-90 
imtermediaie  periodl  (boltoml  instrumentation  at  PAS  The  arrows  show  the  location  of  the  Ps  conversion  discussed  in 
the  icxt  for  both  data  sets.  Note  particularly  the  extreme  amplitude  of  the  Ps  conversion  relative  to  the  first  P  arrival  on 
both  IPN  and  LPN  components 


7 


1942 


Langston:  Scattering  of  Teleteismtc  Body  Waves 


Fig-  7.  Particle  motion  plot  of  the  horizontal  data  for  the  January 
23,  1976,  event.  Data  included  within  the  bracket  are  plotted  below 
with  arrows  occurring  every  0.5  s.  The  Ps  conversion  is  polarized 
almost  perfectly  northward  and  indicates  the  direction  of  dip  of  the 
causative  interface. 


used  to  filter  c(z).  The  filtered  velocity  spectrum  was  inverse- 
transformed  and  scaled  to  a  selected  velocity  variance  and 
mean.  Synthetic  seismograms  were  computed  using  the 
Thompson-Haskell  method  [ Haskell ,  1962], 

Calculations  were  performed  with  a  Gaussian  correlation 
function  as  well,  but  the  exponential  correlation  function 
proved  to  create  more  scattered  arrivals  since  its  spectrum  is 
richer  in  higher  wave  numbers.  A  correlation  length  a  of  1  km 
was  assumed.  For  0.5-Hz  waves  considered  in  this  study,  the 
corresponding  value  of  k,a  is  approximately  equal  to  1,  where 
k,  is  the  vertical  component  of  the  wave  number.  5  wave 
velocities  were  derived  by  assuming  a  Poisson  solid,  and  den¬ 
sity  was  held  constant.  Velocity  parameters  for  the  half-space 
below  the  random  crustal  layer  were  generally  set  to  those  of 
the  lowermost  crustal  layer. 

Two  basic  velocity  models  were  considered  in  the  one¬ 
dimensional  Simula!  ons.  The  first  consisted  of  a  hetero¬ 
geneous  crustal  layer  30  km  thick  with  a  mean  velocity  of  5  5 
km/s  and  velocity  standard  deviation  of  10%.  The  top  pair  of 
wave  forms  in  Figure  10  are  typical  examples  of  the  free  sur¬ 
face  displacements  computed  from  a  number  of  realizations  of 
the  stochastic  parameters.  Figure  9  shows  the  corresponding 
velocity-depth  functions  for  the  lower  and  upper  pair  of  syn¬ 
thetics.  The  incident  wave  time  function  assumed  was  the  time 
derivative  of  the  Gaussian  function  discussed  above  in  the 


source  equalization  section.  A  quick  glance  at  these  synthetics 
ami  the  data  of  Figure  2  show  that  synthetic  coda  levels  are 
significantly  lower  than  those  observed  and  that  the  coda  at¬ 
tenuates  quickly  with  time 

It  might  be  expected  that  a  Moho  with  a  large  velocity 
contrast  would  trap  more  scattered  energy  in  the  crustal  layer 
The  middle  pair  of  synthetics  shows  this  case  for  the  model 
used  to  compute  the  upper  synthetics  but  assuming  a  half¬ 
space  P  wave  velocity  of  8  km/s.  Minor  changes  occur  in  the 
resulting  synthetics.  The  largest  change  is  to  accentuate  the 
Moho  Ps  conversion  by  about  a  factor  of  2  (arrow)  The  coda 
is  largely  unaffected 

Increasing  layer  thickness  tends  to  increase  the  duration  of 
coda.  Increasing  the  velocity  standard  deviation  to  20%  over 
a  layer  60  km  thick  with  a  mean  velocity  of  b  km  s  produces 
synthetic  seismograms  which  start  to  mimic  the  data  'bottom. 
Figure  10).  Large  Ps  conversions  and  reverberations  start  to 
attain  amplitudes  comparable  to  the  direct  radial  P  wave,  and 
coda  duration  superficially  appears  to  agree  with  the  observed 
data. 

Thus  it  appears  that  one-dimensional  velocity  variations  in 
excess  of  20%  over  a  significant  thickness  of  the  lithosphere 
are  needed  to  match  the  receiver  data  at  PAS.  This  would 
suggest  that  a  series  of  layers  which  vary  in  velocity  from 
about  8.4  to  3.6  km/s  from  the  surface  into  the  mantle  over  a 
scale  of  about  I  km  are  present,  which  is  geologically  unrea¬ 
sonable. 

Alternatively,  it  is  possible  that  there  is  significant  body 
wave  to  surface  wa.e  scuitering  in  near-surface  layers  [Damn 
et  al„  1974;  Aki  and  Chouet.  1975;  Dainty  and  Toksoz.  I977] 
This  might  be  expected  on  the  basis  of  the  heterogeneous 
geology  alone.  It  is  also  clear  that  two-  and  three-dimensional 
structures  are  required  in  the  area  to  produce  the  large 
tangential  particle  motions  observed  in  the  data.  Numerical 
experiments  with  two-dimensional  elastic  finite  difference 
methods  show  that,  for  an  equal  amplitude  scattered  field, 
velocity  variations  in  two  dimensions  are  about  half  that 
needed  for  velocity  variations  in  one  dimension  [Mci.auyhtm 
et  dl„  1985] 
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Fig  H  Comparison  of  svnlhctic  radial  wave  forms  for  the  (wo 
models  of  Table  2  with  'he  235  radial  wave  form  slack  The  arrow 
for  the  H-K  (Hadley-kanamort  model)  wave  form  shows  the  location 
of  the  P v  conversion  from  the  Moho  The  arrow  for  ihe  LV/f  (low 
velocity  zone  modell  wave  form  shows  the  location  of  the  Ps  conver¬ 
sion  produced  at  Ihe  base  of  the  crustal  low-velocitv  zone  The  base¬ 
lines  of  the  synthetic  wave  forms  have  been  shifted  for  viewing  pur¬ 
poses 
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Fig  *5.  Velocity  depth  functions  for  'he  Harlley-Kanamori  model 
and  the  LVZ  model  used  in  constructing  the  synthetic  seismograms  of 
Figure  8  Also  shown  are  the  stochastic  one-dimensional  models  used 
in  the  calculation  of  the  synthetics  displayed  in  Figure  10  Note  that 
velocity  and  depth  scales  are  different  among  the  plots. 


Energy  Flux  Models  for  Plane  Wave  Scattering 

Teleseismtc  P  wave  coda  has  been  the  subject  of  a  number 
of  studies  [e.g.,  Aki,  1973;  Dainty  and  Toksoz,  1977;  McLaugh¬ 
lin  et  al„  1985;  Levander  and  Hill.  1985;  Frankel  and  Clayton, 
1986;  Cessaro  and  Butler.  1987].  These  studies  are  closely  re¬ 
lated  to  general  problems  in  coda  generation  due  to  litho¬ 
spheric  structure  in  local  and  regional  data  sets  [e.g.,  Aki. 
1969,  1980;  Aki  and  Chouett,  1975;  Dainty  and  Toksoz.  1977; 
Gao  et  a!..  1983;  Richards  and  Menke,  1983;  Gupta  and  Bland- 
ford.  1983;  Wu  and  Aki,  1985a,  5;  Cessaro  and  Butler,  1987; 
Frankel  and  Wennerberg,  1987;  V idale  and  Helmberger,  1988]. 

A  number  of  techniques  are  available  to  parameterize  the 
coda  level  and  time  decay  based  on  the  Bom  approximation 
of  weak  single  scattering  or  of  diffusion  of  coda  energy  for 
extreme  scattering  [e.g.,  Aki  and  Chouet,  1975;  Aki,  1980; 
Dainty  et  ai.  1974],  Recent  work  has  concentrated  on  simula¬ 
tion  studies  using  finite  difference  acoustic  and  elastic  wave 
propagation  methods  [ Levander  and  Hill,  1985;  Frankel  and 
Clayton.  1986]  which  implicitly  include  the  entire  scattered 
field.  Levander  and  Hill  [1985]  examined  scattering  character¬ 
istics  of  a  rough  boundary  between  a  surface  layer  and  un¬ 
derlying  half-space  and  showed  that  much  of  the  scattered 
held  is  dominated  by  Rayleigh  wave  propagation.  Frankel  and 
Clatton  [1986]  and  Me'  aughlin  e'  al.  [1985]  examined  P-SV 
propagation  in  two-dimensional  random  media  to  examine 
scattering  of  high-frequency  I  f>  1  Hz)  seismic  waves  Subse¬ 


quently,  Frankel  and  Wennerberg  [1987]  developed  a  simple 
theory  based  on  previous  finite  difference  simulations  to  pa¬ 
rameterize  coda  levels,  scattering  attenuation,  and  intrinsic  at¬ 
tenuation  for  two-  and  three-dimensional  scalar  wave  helds. 

The  success  of  a  number  of  receiver  function  studies  in  de¬ 
termining  plane-layered  crustal  and  upper  mantle  structure 
indicates  that  the  scattered  wave  field  may  be  thought  of  being 
composed  of  a  “coherent"  contribution  from  Ps  conversions 
and  reverberations  from  discrete  interfaces  and  a  “stochastic" 
contribution  from  smaller-scale  heterogeneities.  The  coherent 
field  can  be  seen  over  a  large  solid  angle  of  ray  paths  The 
stochastic  field  changes  quickly  with  ray  parameter  and  ray 
azimuth.  Examples  of  the  stochastic  field  are  variation  in 
tangential  P  wave  first  motions  over  the  events  of  the  235 
stack  in  Figure  3 b  as  well  as  coda  arrivals  with  long  lapse 
times  from  the  first  arrival.  The  success  of  any  deterministic 
receiver  function  study  depends  critically  on  the  coherent  field 
being  dominant.  However,  the  incoherent  field,  which  is  usu¬ 
ally  ignored  in  such  studies,  also  contains  statistical  infor¬ 
mation  on  the  degree  of  heterogeneity  in  the  structure  which 
may  be  very  useful. 

A  heuristic  approach  will  be  used  here  to  develop  an  oper¬ 
ational  theory  appropriate  to  the  three-component  receiver 
data.  The  purpose  of  doing  this  is  to  empirically  compare 
coda  levels  and  decay  between  different  receivers  for  classifi¬ 
cation  purposes  and  to  suggest  avenues  of  research  that  will 
address  the  actual  wave  propagation  problems.  This  heuristic 
approaed  will  „,so  be  used  to  quantify  the  differences  oetween 
one-dimensional  structure  coda  development  and  coda  ob¬ 
served  in  the  data. 

A  useful  method  of  parameterizing  the  P  coda  can  be 
derived  following  Frankel  and  Wennerberg  [1987],  They 
examined  scalar  two-dimensional  finite  difference  simulations 
and  suggested  that  scattered  energy  behind  a  cylindrical  or 
spherical  wave  front  distributes  itself  uniformly  over  the 
volume  behind  the  wave  front.  Aki  and  Chouet  [1975]  arrived 
at  the  same  conclusion  when  examining  coda  of  regional 
earthquakes.  This  simple  assumption  yielded  useful  formulae 
for  coda  level  and  decay  in  cases  of  strong  multiple  scattering 
as  well  as  the  limiting  case  of  weak  single  scattering. 

Consider  first  a  source  of  plane  waves  which  radiates  two 
oppositely  propagating  plane  waves  in  a  scattering  whole 


TABLE  2.  Plane-Layered  Crustal  Models 


Layer 

V 

km's 

V„ 

km/s 

Density. 

g/cmJ 

Thickness, 

km 

1 

5.5 

Hadley  and  Kanamori  Model 

3.0  2.6 

4.0 

2 

6.2 

3.5 

2.7 

16.0 

3 

6.8 

3.8 

2.85 

11.0 

4 

7.8 

4.5 

3.1 

10.0 

5 

8.3 

4.8 

3.35 

1 

5.5 

Lon-  Velocity  Zone  Model 

3.0  2.6 

4. 

2 

6.2 

3.5 

2.7 

6. 

3 

6.5 

3.75 

2.7 

5. 

4 

5.6 

7  23 

2.7 

6. 

5 

6.6 

3  81 

2.6 

•> 

6 

7  2 

4  16 

2.9 

7 

7.4 

4.27 
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Fig  10.  Typical  plane  wave  synthetic  seismograms  computed  for  three  realizations  of  one-dimensional  sUK.h.i>iu 
velocity  structure  Vertical  and  radial  displacement  compo  nents  arc  shown.  The  top  traces  are  for  a  mode!  containing  a 
laver  .'0  km  thick  with  an  average  velocity  of  5.5  km  s  and  a  standard  deviation  10";.  of  the  average  The  center  pair  •»* 
synthetics  are  for  the  same  crustal  model  as  the  top  pair  hut  witn  the  addition  of  a  high  velocity  1 8.0  km  >)  half-pace  l  he 
arrow  points  to  the  Moho  Ps  conversion  on  the  radial  component  The  lower  pair  of  synthetics  we  e  computed  i'-tu  .! 
layer  thickness  of  69  km.  average  velocity  of  6  km  s.  and  velocity  standard  deviation  of  20  •  Pn  converMotis  «no 
reverberations  start  to  attain  amplitudes  seen  m  the  PAS  data 


space  The  total  instantaneous  energy  or  power  Er  is  given  by 
the  sum  of  the  direct  wave  power  ED  and  the  coda  power  Ec  . 

Er-t:„+nt. 

Specify  that 


£d  =  £,t 


...  U. 


Ci 


where  £„  is  the  observed  direct  wave  power  and  ••  -  ,  The 
factor  of  2  comes  from  the  fact  that  iw.'  pi. me  waves  are 
propagating  in  the  medium  and  contribute  to  the  scattered 
held  Plane  wave  propagation  theory  is  used  to  obtain  the 
estimate  of  direct  wave  power.  First.  consider  the  integral  of 
me  square  of  the  ground  velocity  -tin 


where  l  is  time,  u>  is  circular  frequency,  and  Q ,  is  the  qualitv 
factor  for  attenuation  due  to  wave  scattering.  Substitution  of 
(2) into  (I)  gives 

Ec  =  fc'^1  -  e  “"u'l  (.11 

Coda  amplitude  Ar  is  related  to  the  coda  power  density 
through  the  principal  assumption  that  the  coda  energy  distrib¬ 
utes  itself  uniformly  behind  the  two  propagating  plane  waves 
First.  we  have 


lD  =  j  -»-(/)  M  ’M 

Times  r,  and  /.  bound  the  direct  wave  arrival  and  are  e>u- 
m.tled  from  ihe  data.  The  direct  wave  power  is  therefore 

i i'En  -  ilr,  PS  i9l 

Substitution  of  |9)  into  ill  and  of  (’!  into  ibl  gives 

■1(.  -  i/n:  •  r1  ’it . :'M  -  -  -  (Id 


1,  - 


(41 


where  d  is  a  scaling  factor  For  P  and  .S  plane  waves  d  =  p  1  3 
where  />  is  density  If  <5S  is  the  unit  plane  wave  area,  r  the 
propagation  distance  from  source  to  receiver,  and  n  P  wave 
velocity,  then  the  coda  power  density  is 


Kj.  =  _I <_  =  T, 

I  2 r  o.S'  2if  o.S 


1 5  > 


where  l  is  volume  Using  (41.  (4).  and  (51  we  obtain 
J{F.,)'  -(l-c  °'ll  ; 

■1(  =  ,  ,  ,  - T-S -  (61 

(2zrl  -  l.'S)  - 

F. ,  can  be  estimated  using  the  observed  dncci  wave  amplitude 
and  correcting  it  for  attenuation  through  the  scattering 
medium  Thus 


l7l 


\s  Frmikil  urn/  UV.i.ierhcv  1  •a']  show,  the  cited  of  at¬ 

tenuation  due  to  scattering  determines  the  uutiai  level  of  the 
coda  scattered  from  the  direct  wave,  nut  coda  decay  with  nme 
is  mainly  controlled  by  the  time-dependent  increase  m  volume 
behind  the  wave  front.  In  this  plane  wave  ease  the  t  ~  ‘  :  de¬ 
pendence  is  due  to  the  linearly  increasing  volume  between  the 
two  oppositely  propagating  plane  waves 

Relation  tl)  is  useful  for  describing  the  coda  for  a  plane 
wave  propagating  in  a  thick  layer  where  coda  lapse  limes 
.wave  arrive  time  relative  to  the  direct  wavel  are  less  than  the 
arrival  time  of  waves  which  interact  with  the  lower  boundary 
of  the  layer,  assuming  the  observation  point  is  at  the  surface 
This  corresponds  to  lapse  times  of  less  than  5  s  for  receivers 
on  typical  continental  crust  A  more  appropriate  model  for 
receiver  function  coda  is  scatlertng  in  a  heterogeneous  layer 
overlying  a  homogeneous,  isotropic  half-space  In  this  situ¬ 
ation  a  vertically  propagating  plane  wave  sweeps  through  the 
layer  once  on  us  way  to  the  receiver,  reflects  from  the  free 
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"face,  and  sweeps  through  yet  another  time  on  its  way  back 
to  the  half-space.  Energy  is  scattered  from  the  plane  wave  into 
coda  energy. 

A  variation  of  this  problem  was  studied  by  Dainiy  et  al. 
[1974]  and  Dainty  and  Toksoz  [1977],  where  they  assumed 
that  scattering  in  the  layer  followed  solutions  to  the  diffusion 
equation.  Also,  assuming  that  all  energy  within  the  layer  was 
scattered  energy,  they  obtained  the  following  analytic  solution 
for  the  scattered  energy  field  at  the  free  surface  m(r)  (assuming 
no  intnnsic  attenuation): 


tn(t) 


1  y  a’  C0S  a"  cx  ( 
h  .f,  la,  4-  sin  2 a,  13  \  4h2  / 


111) 


where  h  is  layer  thickness  and  is  the  vertical  diffusivity  of 
energy  through  the  boundary  of  the  layer  into  the  half-space 
The  coefficients  u.  are  found  as  solutions  of  the  following 
equation 

a,  tan  a,  -  dhvit,,  (12) 

where  v  is  the  seismic  velocity  of  the  half-space. 

The  dominant  term  of  (II)  for  long  lapse  times  and  a  high 
vertical  diffusion  rate  can  be  shown  to  be  for  n  =  1.  Thus  for 
cases  wnere  uirtusion  of  scattered  energy  occurs  quickly,  coda 
energy  decays  like 


m(t)  x  e 


(13) 


where  y  =  c^u, 2  4 h1  This  behavior  can  be  incorporated  into  a 
hybrid  model  containing  aspects  of  plane  wave  propagation  in 
a  layer  with  assumption  of  homogeneity  of  the  scattered  field 
within  the  layer. 

Consider  a  horizontal  scattering  layer  of  thickness  h  overly¬ 
ing  a  homogeneous  and  isotropic  half-space.  A  vertically 
propagating  plane  wave  is  incident  from  below,  passes 
through  the  layer,  reflects  from  the  free  surface,  and  passes 
back  through  the  layer  into  the  half-space.  The  total  power  in 
the  system  can  be  written  as 


Er  =  £„  +  £c  +  £0  (14) 

where  a  new  term  £0  has  been  introduced  to  describe  the 
amount  of  instantaneous  energy  which  diffuses  out  the  bottom 
of  the  layer  at  the  expense  of  the  coda  instantaneous  energy 
Ef  On  the  basis  of  the  behavior  of  ( 11 )  above  I  assume  that 

dEJdt  =  y  Ec  (15) 

which  by  simple  integration  yields 

£0  =  E^e"  —  I )  (16) 

w  here  it  is  specified  that  E0  =  0  at  l  =  0.  As  before,  the  power 
in  the  direct  wave  after  interacting  with  the  layer  is 


E0  =  Erf"2"0, 


(17) 


where  the  factor  of  2  in  the  exponent  comes  from  the  wave 
passing  twice  through  the  layer.  Obviously,  (14)  is  appropriate 
for  short-duration  direct  waves  and  times  greater  than  2tt. 
Substituting  (1 7)  and  ( 16)  into  1 14)  gives 


E,  =  t^l  (18) 


Also,  recognizing  that  the  volume  swept  by  the  plane  wave  is 
now  h  ()S.  the  coda  power  density  becomes 


(19) 


As  before,  the  total  instantaneous  energy  available  to  the 
system  can  be  estimated  from  the  observed  direct  wave  power 
Eg  by 

Er  =  EDe'“‘“Q‘  (20) 

and  using  (9i  gives 

A,  -  (/0'  2  (/  2je‘-'-2tftl  -  e-“2"V  V  "  2  (21) 

Note  that  this  form  for  coda  amplitude  looks  superficially 
the  same  as  that  in  (10)  except  for  the  exponential  factor  of 
time  in  the  numerator  of  (21).  Indeed,  the  time  decay  of  the 
coda  is  controlled  entirely  by  this  factor.  If  y  =  0  is  assumed 
so  that  no  coda  energy  can  diffuse  out  of  the  layer,  then  the 
coda  level  is  constant  for  all  time,  consistent  with  the  plane 
wave  assumption  of  a  packet  of  energy  being  homogeneously 
dispersed  throughout  the  layer.  Thus  the  decay  of  the  coda 
field  is  functionally  equivalent  to  the  leading  term  for  the 
formal  solution  of  111)  and  (13).  particularly  considering  that 
energy  is  proportional  to  the  square  of  amplitude 

Anelastic  attenuation  can  be  included  in  relations  1 10)  and 
|21 )  as  the  factor 

2  *  «■>**  ‘Qip  ~  2c?/ 

where  Q,  is  the  intrinsic  attenuation  of  the  medium  The  first 
exponential  in  this  factor  is  the  correction  factor  to  determine 
total  energy  from  the  direct  wave  and  the  second  exponential 
gives  the  attenuation  of  coda  amplitude.  Note  that  the  effect  of 
coda  energy  diffusing  out  of  the  layer  given  in  (21)  is  exactly 
the  same  as  intrinsic  attenuation.  It  would  be  expected  that  it 
would  be  very  difficult  or  impossible  to  separate  the  two  ef¬ 
fects  uT  practice  using  the  teleseismic  coda  data. 

An  implicit  assumption  in  developing  (10)  or  (21)  is  that  the 
scattered  field  is  of  the  same  wave  type  as  the  primary  These 
equations  are  appropriate  for.  say,  the  scattered  pressure  field 
from  an  incident  P  wave.  Even  for  simple  one-dimensional 
layered  structures,  much  of  the  scattered  field  is  composed  of 
P-to-S  conversions.  For  two-  and  three-dimensional  structures 
there  ir  evidence  that  much  of  the  scattered  field  seen  at  the 
surfac.  is  composed  of  low  group  velocity  surface  waves 
[Dainty  et  al..  1974;  Aki  and  Chouet.  1975;  Levander  and  Hill. 
1985]  Thus  there  is  a  procedural  problem  of  relating  observed 
coda  wave  amplitude  to  energy  since  the  wave  type  contained 
in  the  coda  must  be  known  beforehand.  In  principle,  it  is 
possible  to  directly  infer  the  energy  contained  in  a  wave  field  if 
strain  observations  are  avai1  lble.  However,  three-component 
displacement  data  cannot  be  used  without  assumptions  on 
wave  type 

Recognizing  these  limitations,  equations  (10)  and  (21)  are 
used  as  guides  to  the  analysis  of  the  three-component  data. 
These  equations  will  be  useful  in  parameterizing  relative  levels 
of  coda  and  coda  decay  between  isolated  recievcrs  but  are 
clearly  deficient  in  addressing  all  of  the  scattering  mechanisms 
which  are  probably  important  in  teleseismic  coda  devel¬ 
opment 

Some  operational  aspects  of  examining  coda  decay  are  pat¬ 
terned  after  previous  studies  [e  g.,  Richards  and  Menke.  1983. 
Frankel  and  Wennerherg,  1987],  Observed  three-component 
data  for  a  single  event  aie  fir^t  narrow-band-pass-filtered  with 
a  Butterworth  recursive  filter  in  the  forward  and  backward 
directions  The  two-pole  filter  used  had  corner  frequencies  of 
0  25  and  I  Hz  so  the  following  results  are  appropriate  for  0.5 
Hz  waves  Once  the  data  were  filtered  the  intensity  of  the 
direct  P  wave  at  0?  Hz  was  estimated  by  squaring  the  signal. 
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Fig.  II  Envelope  slack  of  the  ihree-componem  data  of  four 
earthquakes  (see  Table  II  recorded  at  PAS  The  envelope  is  shifted  10 
s  for  display  purposes.  Lines  show  predicted  coda  levels  for  assump¬ 
tion  of  scattering  Q t  of  100,  200.  .TOO.  and  400  using  the  whole  space 
energy  flux  model  (equation  (101).  Scattering  Q ,  at  PAS  is  approxi¬ 
mately  200  to  300,  but  the  coda  appears  lo  decay  slightly  faster  than 
predicted  by  the  model. 


choosing  /,  and  r,  (equation  (8)1  from  the  duration  of  large 
motions  on  the  vertical  component,  and  integrating  over  this 
time  interval.  The  power  of  the  direct  wave  was  estimated 
using  all  three  components  of  motion  over  the  time  interval 
inferred  from  the  vertical  component.  The  integral  of  the 
squared  velocity  used  for  equations  (10)  or  (21)  was  the  square 
root  of  the  sum  of  the  squares  of  the  integrals  from  each 
component  of  ground  motion.  Each  component  was  then 
scaled  by  dividing  the  square  root  of  this  total  squared  veloci¬ 
ty  integral.  The  envelope  of  each  component  was  then  com¬ 
puted  by  forming  the  analytic  signal  [ Farnbach ,  1975]  and 
taking  its  modulus.  The  total  coda  time  series  for  one  event 
was  then  found  by  summing  the  squares  of  the  envelopes  of 
the  three  components  at  each  time  point  and  taking  the 
square  root  of  the  sum.  Resulting  coda  envelopes  for  separate 
events  were  then  averaged  to  obtain  a  better  estimate  of  coda 
level. 

Figure  1 1  shows  the  results  of  this  process  using  four  deep 
earthquakes  recorded  at  PAS  (Table  l|  Deep  events  were 
chosen  to  avoid  contamination  by  near-source  scattering  ef¬ 
fects.  The  observed  levels  of  coda  are  very  high  Indeed,  an 
examination  of  the  raw  data  (e  g..  Figure  21  shows  that  much 
of  the  coda  comes  from  the  horizontal  components.  Theoreti¬ 
cal  curves  computed  using  equation  (101  arc  superimposed  on 
the  coda  decay  curve  in  Figure  1 1  and  show  that  an  apparent 
scattering  Qt  of  200  to  300  is  required.  The  coda  time  decay 
appears  to  be  very  slow  and  is  roughly  consistent  with  r~‘  2 
found  in  this  model  of  scattering 

Figure  12  demonstrates,  however,  that  the  simple  one- 
dimcnsional  simulations  are  not  consistent  with  coda  decay 
following  equation  ( 1 0).  The  coda  curve  for  the  "20''.."  model 
was  constructed  by  stacking  10  vertical  and  horizontal  com¬ 
ponent  realizations  of  models  which  had  a  velocity  standard 
deviation  of  20"-..  and  a  layer  thickness  of  hO  km.  TK'  “10%“ 
curve  was  obtained  by  stacking  nine  vertical  and  horizontal 
component  realizations  for  models  which  had  a  velocity  stan¬ 
dard  deviation  of  10%  and  a  layer  thickness  of  30  km  The 
observed  coda  decay  is  linear  on  the  logarithm  plot  and  falls 
off  much  faster  than  implied  by  1 1 01  The  linear  lall-off  is 
consistent  with  the  scattering  layer-over  half-space  model 


where  coda  energy  diffuses  out  of  the  layer  into  the  half-space 
governed  by  equation  (15).  The  one  dimensional  simulations 
included  no  effect  of  anelastic  attenuation 

Figure  12  also  shows  least  squares  linear  his  to  the  coda 
decay  to  obtain  Q,  and  y  in  equation  (21)  The  slope  of  the 
log-coda  curve  yields  /,  and  the  zero  time  intercept  can  bo¬ 
used  to  solve  directly  for  Qr  The  standard  deviation  of  the 
least  squares  fit  was  also  used  to  estimate  allowable  Q ,  vari¬ 
ation  by  adding  and  subtracting  the  >tandard  deviation  from 
the  zero  intercept  time  to  find  a  lower  and  upper  bound  of  Qt, 
respectively.  These  values  are  displayed  in  Table  3.  Coda  from 
the  simulations  show  that  the  diffusing  layer  model  correctly 
predicts  the  form  of  coda  decay  although  the  model  does  not 
formally  treat  the  scattering  mechanism  of  P  to  S  conversions 
within  the  layer.  The  decay  rate  is  very  sensitive  to  the  velocity 
standard  deviation,  but  Qt  estimates  are  surprisingly  the  same, 
within  the  error  of  determination. 

It  is  interesting  to  compare  results  for  PAS  with  those  from 
another  station  to  get  an  appreciation  for  the  level  of  scatter¬ 
ing  implied  oy  me  data.  Tnree  deep  evem>  recorded  on  the 
broadband  DWWSSN  system  at  State  College,  Pennsylvania 
(SCP)  were  analyzed  in  the  same  wav  Event  parameters  can 
be  found  in  Table  4.  and  the  data  are  displayed  in  Figure  13 
The  Benioff  1-90  and  intermediate  period  DWWSSN  systems 
are  sufficiently  similar  for  the  purposes  of  this  comparison, 
particularly  since  the  same  band-pass  filter  was  used  on  the 
data. 

Figure  14  compares  the  coda  decay  curve,  for  PAS  and 
SCP.  Structure  under  SCP  is  seen  to  be  simple;  than  that  at 
PAS  [ Langston  and  Isaacs ,  1981,  Ammon,  personal  communi¬ 
cation,  1988]  and  gives  rise  to  lower  amplitude  Ps  conversions 
as  well  as  coda  Coda  decay  for  SCP  is  twice  as  fast  as  that 
observed  for  PAS  (Table  3).  Q,  is  found  to  be  lower  for  PAS 
with  use  of  equation  f21).  giving  a  value  of  239  compared  to 
582  for  SCP. 
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Kig  12  C  oda  envelope  stacks  for  the  > >ne-dimensional  simula¬ 
tion  The  "10  "  envelope  is  the  stack  of  1  >  vertical  and  radial  svn- 

t  he  tic  seismograms  produced  by  nine  realizations  of  the  ,'0-km-ihick 
layer  model  with  average  velocity  of  *5  kin  >  and  velocity  standard 
deviation  of  10  The  "20  envelope  is  the  stack  of  20  vertical  and 
rjdul  synthetic  seismograms  produced  hv  !0  realizations  of  the  t>0 
km  thick  layer  model  with  average  velocitv  of  6  km  s  and  velocity 
•standard  deviation  of  20  The  straight  hnes  are  least  squares  hts  of 
the  coda  values  are  those  inferred  from  the  zero  lapse  time  mtcr- 
ewpt  of  the  hne.tr  tits  (see  (evu  Note  :ha*  viaJj  decay  in  the  one- 
dimetisional  simulations  agree  with  the  assumption  that  soda  energy 
follows  a  diffusion  law  for  leaking  into  the  h.iil-'pase  (  i*da  decay  »> 
controlled  entirely  hv  this  process  in  the  one -‘Jimen  signal  simulations 
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TABLE  3.  Q,  and  y  Determinations  lor  Coda  Stacks 


Coda  Stack 

Least  Squares  Fit 

Intercept  Slope 

Standard 

Deviation 

V, 

y 

High  Q 

Low  Q 

10^-  simulation 

-0.970 

-0.050 

0.079 

544 

0.229 

7K4 

380 

20^  simulation 

-0  985 

-0  0098 

0  057 

-'S4 

0.045 

'61 

450 

PAS  data 

-0  789 

-0.0047 

0.075 

239 

0.022 

337 

169 

SCP  data 

-0.984 

-0.0092 

0  105 

5X2 

0.043 

944 

359 

Scattering  Q  determinations  were  found  using  equation  (21)  for  a  )l)-km-lhick  scattering  laser  with 
an  average  P  velocity  of  6  km/s. 


Discussion 

The  scattering  layer-over  half-space  model  reproduces  the 
principal  behavior  of  the  one-dimensional  st'ucture  simula¬ 
tions  (Figure  12)  and  is  consistent  with  coda  decay  in  the  PAS 
and  SCP  data  The  simple  assumptions  of  homogeneity  of  the 
coda  field  and  diffusion  of  energy  into  the  half-space  seem  to 
describe  the  basic  mechanisms  of  coda  formation  and  is  con¬ 
sistent  with  previous  observations  of  the  behavior  of  data  and 
two-dimensional  simulation  studies 

4ki  and  Chouet  [1975]  estimated  the  diffusivity  of  the  litho¬ 
sphere  in  Japan  and  California  using  a  diffusion  model  of 
coda  formation  applied  lo  local  earthquake  data  They  found 
high  diffusion  rates  having  the  effect  of  homogenizing  the  coda 
field  behind  the  wave  front.  Frankel  and  Wennerheri/  [1987] 
took  these  ideas  further  by  examining  the  coda  field  in  finite 
difference  simulations  and  constructing  a  simple  energy  flux 
theory  to  explain  the  formation  of  coda  Although  the  as¬ 
sumptions  of  homogeneous  coda  and  diffusive  energy  flow 
across  the  layer  boundary  are  reasonable,  the  actual  mecha¬ 
nisms  of  coda  formation  are  not  directly  addressed  in  an  equa¬ 
tion  like  (21),  which  leads  to  the  problem  of  estimating  coda 
energy  from  an  unknown  wave  field. 

Much  of  the  coda  in  the  one-dimensional  simulations  is  a 
product  of  P  to  S  conversions  and  reverberations.  The  energy 
scattered  into  S  waves  is  obviously  a  function  of  ray  parame¬ 
ter  As  the  direct  wave  incidence  angle  increases,  more  P  to  S 
conversions  will  occur  This  can  be  verified  directly  by  calcula¬ 
tion  but  can  be  seen  in  the  behavior  of  the  conversion  coef¬ 
ficient  at  a  boundary  Thus  it  can  be  expected  that  coda  fall- 
off  and  levels  will  change  for  waves  of  different  incidence  angle 
if  one-dimensional  structure  is  appropriate  P  to  S  scattering 
in  two-dimensional  structures  is  more  complex  [Frankel  and 
Clayton,  1986;  McLauqhin  el  at..  1985]  but  appears  to  become 
less  sensitive  to  incidence  angle  P  and  S  to  Rayleigh  scatter¬ 
ing  is  piuoubly  a  major  component  of  the  coda  field  at  rela¬ 
tively  low  frequencies  (  <  1  Hz)  [,4ki  and  Chouet.  1975]  These 
scattering  mechanisms  may  control  the  coda  formation  in  the 
data  presented  here.  In  terms  of  the  application  of  equation 
|21)  the  problem  amounts  to  estimating  d,  the  scaling  factor 


relating  energy  density  to  wave  amplitude  in  equation  14). 

Even  for  one -dimensional  structure,  d  depends  on  incidence  j 

angle  and  includes  the  free  surface  receiver  functions  [eg,  ! 

Helmberqer,  1468]  It  is  ol  some  interest  to  examine  ihe  energy 
partitioning  in  ihe  coda  of  the  simulations  and,  making  some 
simple  assumptions,  the  partitioning  seen  in  the  data 

As  an  approximation,  consider  ihe  coda  power  being  com-  j 
posed  of  S  wave  £(.  and  P  w  ave  Ec  powers 

£c  =  £c  *  ECf  1 22)  j 

Also  deline  the  energy  partuoning  coefficient  by 

r  =  ^  1 23) 

f-c. 

For  plane  wave  propagation 

t'c.  =  Pffl, 

1 24) 

*•«’.  = 

where  /,  and  I p  are  the  estimated  integrals  of  squared  velocity 
for  S  wave  and  P  wave  motions,  respectively  The  5  wave 
velocity  is  given  by  /i  For  one-dimensional  structure  models 
and  for  incident  P  waves  of  small  incidence  angle.  S  waves 
occur  primarily  on  the  horizontal  component  and  P  on  the 
vertical.  The  respective  wave  integrals  can  therefore  be  directly 
estimated  using  (8)  by  performing  the  integration  over  the 
filtered  and  squared  wave  form  from  the  end  of  the  direct 
wave  arrival  to  some  reference  ume  in  the  coda.  This  was 
done  for  the  wave  forms  obtained  from  the  1 0° ..  and  20"., 
one-dimensional  simulation  models  Both  models  give  similar 
results  where,  lor  the  ray  parameter  considered  (0  06  s  km)  in 
constructing  the  synthetics,  roughly  70°.,  of  the  scattered 
energy  occurs  as  .S'  wave  energy.  I.A  value  of  0.7  +  0.3  was 
obtained  for  both  simulations  using  the  individual  wave  forms 
of  each  model  realization  I  Changing  the  ray  parameter  to  0  04 
s  km.  appropriate  for  source  distances  of  85  .  reduces  the  S 
wave  energy  to  50  and  less  for  the  synthetics.  The  free  sur¬ 
face  effect  is  assumed  to  be  the  same  for  both  P  and  S  for 
these  low  angles  of  incidence 


TABLE  4  Event  Parameters  for  State  College  Data 


Origin 


Date 

Time, 

I'T 

Latitude. 

deg 

Longitude. 

deg 

V/„ 

Depth. 

km 

Distance. 

sL'g 

Back  A/imuth. 
deg 

Dec  21.  1987 

1214  18 

28  2S 

M  :w 

6*02 

-o  0 

1^3 

Aug  8.  1985 

1  #>3^-58 

6  IS 

t  M  4f 

<  - 

144  0 

Ml 

Aug  12.  198  s 

<)4Vi;4? 

t  os 

IP  .21. 

41 

|4l  ,H 

U4 

194* 
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Fig.  13  Three-component  data  recorded  at  SCP  on  the  DWWSSN  intermedtate-penod  system  iTable  4i  Long-period 
noise  seen  on  some  horizontal  wave  forms  teg.  August  12.  1985)  was  largely  removed  by  the  band-pass  biter  used  in  this 
study 


Small  differences  in  correcting  the  data  for  wave  propaga¬ 
tion  effects  are  of  little  consequence  to  this  discussion  since 
wave  types  in  the  observed  coda  data  are  largely  unknown. 
We  treat  the  observed  data  in  the  same  way  where  the  5  and 
P  wave  integrals  are  defined  as 


where  the  subscripts  Z,  V.  and  E  denote  the  component  of 
ground  motion  Assuming  only  5  and  P  wave  partitioning  in 
the  PAS  data  yields  a  partitioning  coefficient  of  1.7  <-0.4.  a 
factor  of  2  3  greater  than  expected  compared  to  the  one- 
dimensional  simulations  This  result  is  consistent  with  the 
coda  being  composed  of  low  group  velocity  surface  waves 
scattered  from  incident  P  and  5  waves  Instantaneous  energy 
will  be  proportional  to  the  group  velocity  so  that  assuming  a 
higher  S  velocity  in  (25)  will  cause  the  energy  to  be  overesti¬ 
mated  These  observations  are  consistent  with  observations  of 
the  coda  at  arrays  [,4Ai  1973.  -Hi  arul  Chouel.  1975]  and  from 
theoretical  wave  propagation  calculations  fe.g..  Levander  and 
Hill,  1 985(]  Powell  and  Meltzer  |_|9S4]  also  found  direct  evi¬ 
dence  for  a  high  level  of  scattering  under  southern  California 
in  their  study  of  coherency  across  the  C  altech-C  S  Geological 
Survey  large  seismic  array  (SCARLET I. 

It  thus  appears  that  coda  level  and  coda  decay  at  PAS  is 
inconsistent  with  plausible  one  dimensional  Earth  models 
The  observed  data  show  slow  coda  decay,  implying  a  rela¬ 
tively  long  dwell  time  of  coda  energy  in  the  crust  as  well  as 
high  amplitudes  The  high  amplitudes  are  consistent  with  the 
coda  Held  being  primarily  composed  of  scattered  surface 
waves  Even  the  coda  decay  seen  at  SCP  implies  unreasonable 
one-dimensional  structure,  since  the  data  imply  virtually  the 
same  attenuation  and  decay  as  the  20  velocity  model  iTable 
3) 

These  aspects  of  the  receiver  function  data  can  be  routinely 
quantified  in  other  data  sets  to  motivate  an  interpretation  of 
structure  under  a  receiver  If  the  receiver  function  coda  data 
show  tendencies  that  are  inconsistent  with  simple  one- 
dimensional  models,  then  inversion  of  selected  phases  at  long 
lapse  times  1  >  10  si  from  the  first  arrival  or  inversion  of  the 
entire  wave  form  becomes  suspect  I  his  result  is  not  surprising 
when  one  simply  looks  at  the  anomalies  contained  within  the 


data,  but  the  theoretical  trealment  presented  here  can  help 
quantify  both  the  gross  characteristics  of  the  observed  data 
and  the  justification  of  a  particular  modeling  strategy 

Ps-P  arrival  times  suggest  that  a  major  discontinuity  occurs 
in  the  midcrust  under  PAS  Although  the  polarity  of  radial 
and  tangential  Ps  is  consistent  with  the  interface  dipping 
northward  under  the  San  Gabriel  Mountains,  dip  is  of  the 
order  of  10  or  less,  and  the  S  wave  velocity  contrast  must  be 
unusually  large  I  I  km  s)  Qualitatively,  a  major  crustal  low 
velocity  zone  can  explain  these  observations,  but  the  ex¬ 
tremely  large  Pv  u  amplitude  ratios  probably  imply  that  other 
factors  are  affecting  the  wave  form  such  as  ray  focusing  [Lee. 
1983 .  Lee  and  Langston,  1983u.b] 

The  northward  dip  of  the  interface  suggests  that  it  is  a 
major  structure  associated  with  the  southward  overthrusting 
of  the  San  Gabriel  Mountains  If  the  amplitude  of  the  conver¬ 
sion  is  due  to  large  velocity  contrasts,  then  a  low  velocity  zone 
is  required  at  midcrustal  to  ueep  crustal  levels.  It  is  interesting 
to  note  that  this  low  velocity  zone  occurs  just  under  the  sets- 
mogente  zone  of  the  region  and  may  he  the  seismic  signature 
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of  the  decoupling  zone  of  upper  crustal  and  lower  crustal 
upper  mantle  microplates  [ Humphreys .  1984;  Wehband  Kana- 
mon,  1985]. 

The  high  velocity  mantle  anomaly  under  the  Transverse 
Ranges  inferred  from  past  studies  may  have  some  effect  in 
perturbing  the  incident  F  wave  ray  path.  Rays  will  tend  to  be 
bent  toward  the  vertical  and  toward  the  interior  of  the  anoma¬ 
ly.  However,  the  expected  perturbation  will  only  be  a  few 
degrees  The  observed  scattering  effects  induced  by  crustal  het¬ 
erogeneity  dominate  the  character  of  the  horizontal  compo¬ 
nents  of  the  teleseismic  wave  forms 

Conclusions 

The  receiver  function  data  set  for  PAS  suggests  that  the 
scattering  is  occurring  in  a  highly  heterogeneous  crust  Broad¬ 
band  Benioff  1-90  data  from  teleseisms  show  anomalous 
tangential  particle  motions  and  a  high-amplitude  coda  which 
decays  slowly  Initial  portions  of  the  radial  and  tangential 
receiver  functions  show  a  coherent  inferred  Ps  conversion 
which  displays  a  polarization  anomaly  of  45  for  most  data 
Using  the  amplitude,  polarity,  and  timing  of  this  phase  seen  in 
stacks  of  the  data  and  from  a  direct  observation  in  an  incident 
PKIKP  phase,  a  high  S  wave  velocity  contrast  (>!  km  si 
interface  is  inferred  at  approximately  20  km  depth.  The  inter¬ 
face  dips  less  than  10  to  the  north  and  appears  to  be  a  major 
structure  associated  with  southward  overthrusting  of  the  San 
Gabriel  Mountains 

Observed  coda  level  and  decay  was  examined  using  two 
methods  One  was  direct  simulation  of  one-dimensional  sto¬ 
chastic  structures.  Plane  wave  synthetic  seismograms  were 
computed  for  random  plane-layered  models  with  an  ex 
ponential  correlation  function  and  with  10%  and  20%  stan¬ 
dard  deviations  in  velocity  The  PAS  data  showed  larger  scat¬ 
tering  effects  than  the  simulations,  indicating  that  geologically 
unreasonable  one-dimensional  models  are  required  to  explain 
the  coda  data  The  one-dimensional  models  also  are  obviously 
deficient  in  explaining  the  degree  of  ofT-azimuth  scattering 
seen  in  the  data. 

The  other  method  consisted  of  examining  coda  behavior 
using  an  energy  flux  model  developed  for  a  scalar  plane  wave 
incident  on  a  scattering  layer  over  a  homogeneous  half-space 
A  scattering  layer  model  was  considered  since  it  is  likely  that 
major  velocity  perturbations  are  largely  confined  to  the  crust 
Two  fundamental  assumptions  were  made  to  develop  the 
model  and  were  based  on  previous  empirical  observations  of 
the  behavior  of  earthquake  coda  and  numerical  experiments. 
It  was  assumed  that  the  coda  held  distributes  itself  homoge¬ 
neously  within  the  layer  and  that  coda  energy  diffuses  across 
the  layer-half-space  boundary.  Coda  decay  is  seen  to  be  con¬ 
trolled  entirely  by  the  diffusion  constant  of  cnere.  f.c/.v  across 
the  layer  boundary.  Synthetic  seismograms  from  the  onc- 
dimensional  simulations  show  that  the  simple  energy  flux 
model  explains  the  form  of  coda  decay  One  implication  of 
this  model  is  that  the  diffusion  effect  is  indistinguishable  from 
anelastic  attenuation  Thus  it  is  likely  that  teleseismic  coda 
data  cannot  be  used  to  estimate  local  anelastic  attenuation 

PAS  and  SCP  data  from  selected  deep  earthquakes  were 
unulv/ed  using  this  model  and  it  was  found  that  PAS  had  a 
i  -wet  vaiiering  (J,  l  -  259)  compared  to  SCP  (  -  582)  and  that 
Mu-  vi >da  decay  for  SCP  was  twice  as  fast  as  that  for  PAS  The 
ubiolule  values  of  scattering  f),  obtained  with  t he  model  arc 
-uhic-ci  to  assumptions  on  the  fvpes  of  waves  contained  with 
rtie  wave  held  arid  piotiuhlv  repiesein  lower  hounds  In  die 


actual  Q  values  The  comparison  between  the  two  stations 
shows  i hat  scattering  is  lower  in  a  tectonically  quiescent  area 
with  less  variable  geology  as  expected. 

An  analysis  of  energy  in  the  horizontal  and  vertical  compo¬ 
nents  of  ihe  one-dimensional  synthetics  and  the  PAS  data 
suggests  that  much  of  the  energy  contained  in  the  observed 
coda  is  from  scattered  surface  waves. 


Appendix 

Nominal  instrument  constants  for  the  Bemoff  1-90  system 
arc  pendulum  period  of  I  s,  galvonometer  period  of  90  s. 
damping  constants  of  I  ftu  both  the  seismometer  and  gal¬ 
vonometer.  a  coupling  constant  of  0.05.  and  magnification  of 
5000  iH.  K.uumori.  personal  communication.  1987i  A  cali¬ 
bration  of  ihe  system  was  started  in  1962  but  was  never  totally 
completed  Calibration  of  the  vertical  component  showed  a 
peak  magnification  ol  2700.  10'.,  under  nominal  specifi¬ 
cations  Experience  wiih  ihe  system  suggests  that  instrument 
constants  arc  good  to  about  30".,.  Because  calibration  of  the 
instruments  can  affect  the  results  of  rotation  of  the  data  and 
the  source  equalization,  it  is  of  some  interest  to  examine  the 
results  of  errors  in  the  instrument  constants. 

Assuming  that  the  receiver  response  is  ideal  and  consists  of 
motions  confined  to  the  sagittai  plane  containing  the  ray. 
vector  rotation  of  the  horizontal  displacement  components  to 
obtain  radial  and  tangential  ground  motions  yields 


!i„i n  =  Ru\  •  [ijrl  cos:  9  r  utn  sin:  d] 
’irUI  -  Pin  •  [/.in  -  !„(/)]  sin  20  2 


(All 


where  the  subscripts  R  and  /  denote  radial  and  tangential 
motions,  respectively  ,  and  the  subscripts  n  and  e  denote  north- 
south  and  east-west  motions,  respectively  The  it  is  the  back 
azimuth  angle  to  the  source  from  the  receiver  The  respective 
instrument  responses  are  given  by  i,(il  and  i,(M.  and  Kill  is  a 
common  radial  response  for  plane  layered  structure 

The  largest  errors  in  rotation  occur  when  0  =  45  Tangen¬ 
tial  motion,  in  this  case,  is  caused  by  differences  in  the  instru¬ 
ment  response  of  the  two  components  Clearly,  if  small  differ¬ 
ences  occur  in  ihe  response,  then  the  radial  motions  will  be 
little  affected,  since  the  net  response  will  be  the  average  of  the 
two.  Magnification  errors  will  primarily  show  up  in  iht 
tangential  component  In  the  worst  case  considered  h.re. 
langcnti.il  motions  will  be  approximately  50'..  oi  the  radial 
component,  if  magnification  is  only  known  to  50  .,  and 
should  look  identical  (o  [he  radial  component  in  wave  shape 
The  data  for  P  AS  l Figure  2)  show  extreme  differences  in  wave 
shape  between  the  horizontal  components,  which  cannot  be 
due  to  magnifi-ut.on  errors 

The  equations  for  electromagnetic  seismographs  [HatfiHnra. 
1 95S ]  were  used  to  estimate  the  difference  in  instrument  re¬ 
sponses  it  fit  variations  in  the  instrument  constants  are  as¬ 
sumed  Figure  A  I  displays  tile  results  for  amplitude  spectra 
Changes  of  9)  m  ihe  galvonometer  period  and  damping  and 
ihe  seismometer  period  and  damping  were  assumed  relative  to 
the  nominal  response  I  heoretical  responses  were  calculated  in 
the  frequenev  domain  and  inverse  Fourier-translormed  to 
obtain  impulse  -espouses  Using  equation  (All  as  a  guide,  ihe 
perturbed  responses  were  then  subtracted  from  the  nominal 
response  and  then  F.uirier-iransformcd  to  obtain  the  ttnph 
lu.le  spectrum  I  bus  On-  i.uir  curves  below  dir  nominal  1  '*0 
r.  J.on-.e  -,.-n  *.  1  i'll'.-  \l  ale  tin-  amplitude  specie!  ol  1  tie- 
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Fig  AI.  Amplitude  spectra  of  ihe  nominal  Benioff  1-90  response 
(topi  and  response  differences  (lower  curvesl  assuming  30";,  variation 
in  damping  and  free  periods  of  the  seismometer  and  gafvonome(er  of 
the  system.  Parameters  hs.  hg,  Ts.  and  Tg  are  (he  seismometer  damp¬ 
ing.  galvonometer  damping,  seismometer  period  (in  secondsl.  and  gal- 
vonometer  period,  respectively  See  text  for  explanation 


differenced  impulse  responses  They  can  be  considered  numeri¬ 
cal  derivatives  of  the  instrument  response  if  divided  by  0.3. 

A  change  in  the  galvonometer  period  or  damping  results  in 
a  response  2-3  orders  of  magnitude  lower  than  the  nominal 
response  in  the  band,  centered  about  1  Ha  I  Figure  Al).  Thus  it 
is  not  likely  that  errors  in  these  parameters  will  be  of  any 
consequence  in  the  data.  Changes  of  30%  in  seismometer 
period  produces  a  tangential  impulse  response  about  a  factor 
of  4  lower  than  the  nominal  response  and  looks  nearly  identi¬ 
cal  to  the  nominal  response.  A  change  in  seismometer  period 
appears  as  a  change  in  magnification.  The  tangential  wave 
form  would  differ  by  only  a  constant  compared  to  the  radial 
wave  form.  A  change  in  seismometer  damping,  however,  has 
the  greatest  change  in  the  shape  of  the  spectrum.  A  90  phase 
shift  is  evident  at  1  Hz  which,  for  band-limited  data,  would 
make  the  tangential  motion  appear  Hilbert-transformed  com¬ 
pared  to  radial  motion  Fortunately,  this  response  is  1-2 
orders  of  magnitude  lower  than  the  nominal  response  In  sum¬ 
mary,  plausible  changes  in  the  instrument  constants  for  the 
1-90  system  cannot  explain  the  anomalous  particle  motions 
seen  in  the  data. 

Magnification  errors  cannot  be  discounted.  However,  an 
empirical  test  was  made  by  comparing  the  ratio  of  north- 
south,  east-west,  and  vertical  amplitudes  of  the  first  P  pulse 
observed  in  the  1-90  data  with  that  seen  in  the  30-90  data.  P 
wave  data  for  the  December  28.  1973,  March  23.  1974.  and 
November  29,  1974.  events  were  used  Considering  that  the 
passbands  of  the  two  instruments  are  different,  amplitude 
ratios  of  the  different  ground  motion  components  between 
instruments  were  within  20%  of  each  other 

Finally,  the  data  can  be  used  in  the  test  proposed  by  Lang¬ 
ston  [1979]  to  demonstrate  the  major  off-azimuth  arrivals 
occur  on  both  horizontal  components  for  events  with  different 
back  azimuths.  Figure  A2  shows  the  first  11)  s  of  Ihe  P  wave 
forms  for  the  February  I.  1973  and  December  28.  1973.  events 
(also  see  Figure  2).  Polarities  ha>  e  been  adjusted  to  make  the 
wave  form  comparison  clearer,  and  the  P  waves  have  been 
aligned  in  time  The  arrow  in  the  middle  pair  of  plots  shows 
the  location  of  the  large  Ps  conversion  studied  in  the  main 
body  of  the  paper  It  clearly  occurs  on  ihe  north-south  com- 
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Fig.  A2.  Comparison  of  Benioff  1*90  three-component  data  from 
the  February  1.  1973.  event  (top)  and  the  December  28.  1973.  event 
(bottom!.  Polarities  have  been  reversed  f  -  the  vertical  (Z)  and  east- 
west  iE)  components  of  the  December  28.  1973.  event  for  comparison 
purposes.  The  vertical  components  are  simple,  showing  a  single  im¬ 
pulsive  P  wave.  The  arrow  shows  the  location  of  the  major  Ps  conver¬ 
sion  considered  in  the  this  study  It  occurs  primarily  on  the  N  compo¬ 
nent  for  the  February  1.  1973.  event  and  on  the  E  component  for  the 
December  28.  1973.  event.  Likewise,  it  is  not  obvious  on  the  other 
respective  horizontal  components  (right  side),  showing  that  both  hori¬ 
zontal  instruments  respond  similarly  to  the  same  wave  propagation 
effect 


ponent  for  the  February  I.  1973.  event  and  on  the  east-west 
component  of  the  December  28,  1973.  event  Likewise,  the 
corresponding  east-west  and  north-south  components  (right 
pair  of  wave  forms)  show  similar  wave  forms  between  events 
without  the  major  arrival.  This  comparison  shows  that  both 
horizontal  instruments  behaved  in  a  similar  fashion  for  the 
same  wave  propagation  effect. 
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Abstract 

The  codas  of  long-period  Rayleigh  waves  recorded  at  WWSSN  and 
Canadian  network  stations  in  Western  North  America  from  eight 
underground  explosions  at  NTS  are  examined  in  an  effort  to  sepa¬ 
rate  scattering  and  anelastic  attenuation  effects.  Coda  behavior 
of  0.1  and  0.2  hz  Rayleigh  waves  follows  coda  characteristics 
seen  in  studies  of  short-period  S  waves.  Coda  decay  rate  is  seen 
to  be  a  stable  observation  over  most  stations  in  Western  North 
America  and  is  consistent  with  the  hypothesis  that  backscattered 
surface  waves  from  heterogeneities  contained  within  the  western 
half  of  the  continent  form  the  Rayleigh  wave  coda.  The  basic 
data  observables  of  coda  level  and  decay  are  interpreted  using 
several  plausible  models.  The  single  scattering  model  yields  a 
coda  Q  consistent  with  previously  determined  Rayleigh  anelastic 
attenuation  coefficients.  Separation  of  anelastic  and  scattering 
Q  is  possible  using  an  energy  flux  model  and  shows  that  scatter¬ 
ing  Q  is  one  to  two  orders  of  magnitude  higher  than  anelastic  Q. 
However,  an  energy  flux  model  which  incorporates  a  layer  of  scat- 
terers  over  a  homogeneous  halfspace  shows  that  all  Rayleigh  wave 
attenuation  can  £>e  explained  purely  by  scattering  effects  which 
include  Rayleigh  to  body  wave  conversion.  Coda  can  be  fit 
equally  well  by  these  mutually  incompatible  models.  It  is  not 
likely  that  the  mechanisms  of  scattering  or  anelastic  attenuation 
can  be  addressed  by  coda  observations  of  a  single  homogeneous 


data  set. 
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Introductlon 

Recent  interest  in  characterizing  lateral  heterogeneity  in 
the  lithosphere  has  led  to  increasing  study  of  coda  waves  which 
are  a  direct  result  of  this  heterogeneity.  The  working  model 
developed  by  Aki( 1969 ; 1973 ; 1980a ; 1980b)  and  Aki  and  Chouet(1975) 
that  the  S  wave  coda  consists  of  backscattered  S  waves  and  is 
characterized  by  a  single  scattering  Bom  approximation  has  found 
wide  use  in  data  interpretation  of  local  seismograms  (see  the 
review  article  by  Herraiz  and  Espinosa,  1987).  This  simple 
model  explains  much  of  the  character  of  observed  coda.  What  is 
not  so  evident,  however,  and  is  a  point  of  some  controversy,  is 
the  interpretation  one  makes  of  the  coda  Q,  Qc ,  parameter. 
Strictly  speaking,  Qc  is  the  attenuation  effect  due  to  scattering 
of  the  elastic  wave  by  elastic  heterogeneity.  Empirically,  Aki 
and  Chouet(1975)  and  Aki(1980b)  have  suggested  that  Qc  is  more 
closely  related  to  anelastic  Q.  This  interpretation  was  quanti¬ 
fied  by  Frankel  and  Wennerberg( 1987 )  and  Wu(1985)  who  developed 
energy  flux  models  which  include  anelasticity  and  multiple  scat¬ 
tering  of  scalar  waves.  In  Frankel  and  Wennerberg's  model,  scat¬ 
tering  is  manifest  principally  in  coda  amplitude  levels  while 
anelasticity  dominates  the  coda  decay  with  time. 

In  this  paper,  a  coda  data  set  consisting  of  long-period 
Rayleigh  waves  recorded  within  western  North  America  is  examined 
using  techniques  normally  applied  to  high-frequency  S  wave  codas. 
A  primary  motivation  for  this  unusual  approach  was  to  examine  the 
stability  and  empirical  characteristics  of  coda  amplitude  decay 
of  the  fundamental  mode  Rayleigh  wave  for  possible  use  in  source 
magnitude  and  yield  estimates. 


JO 
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Rayleigh  wave  periods  of  20,  10  and  5  seconds  are  examined 
over  coda  lapse  times  of  1000  seconds  (from  source  origin  time). 
At  these  long  lapse  times  and  for  Rayleigh  wave  group  velocities 
of  approximately  2.5  km/sec,  scattered  waves  comprising  the  coda 
can  originate  from  heterogeneities  distributed  over  western  North 
America  and  the  eastern  Pacific  Ocean.  This  data  set  yields 
information  on  large  scale  heterogeneity  in  the  continental 
crust . 

Rayleigh  coda  will  be  examined  using  three  different  coda 
theories.  The  standard  single  -  scattering  model  of  Aki(1969)  will 
be  used  to  obtain  Qc  and  to  compare  it  with  previous  determina¬ 
tions  of  anelastic  Q  made  from  spectral  amplitude  decay  measure¬ 
ments  of  fundamental  mode  Rayleigh  waves  in  western  North  America 
(Mitchell,  1975).  The  separation  of  anelastic  Q  and  scattering  Q 
will  be  attempted  using  the  energy  flux  model  of  Frankel  and  Wen- 
nerberg(1987) .  Their  model  incorporates  the  effects  of  multiple 
scattering  and  conservation  of  wave  energy  and  may  represent  a 
more  realisti^  wave  propagation  situation  compared  to  the  single¬ 
scattering  model.  A  third  model  based  on  the  energy  flux  model 
but  incorporating  the  effect  of  surface  wave-to-body  wave  conver¬ 
sions  from  scattering  in  a  heterogeneous  layer  over  a  homogeneous 
halfspace  is  also  used  to  determine  if  this  plausible  scattering 
mechanism  can  explain  the  data.  This  last  model  was  developed  by 
Langston(1988)  to  explain  P  wave  coda  generated  by  teleseismic  P 
waves  scattering  under  a  receiver.  It  was  found  that  the  three- 
component  coda  of  teleseismic  receiver  functions  recorded  at 
Pasadena,  California,  had  characteristics  not  explainable  by  the 
single  scattering  model  and  was  also  consistent  with  body  wave- 
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to-surface  wave  conversion.  The  assertion  is  made  that  the  reci¬ 
procal  problem  of  surface  wave-to-body  wave  conversion  is  equally 
important  in  the  scattering  of  fundamental  mode  surface  waves. 

The  results  of  this  study  show  that  it  is  possible  to  fit  the 
coda  data  equally  well  with  each  theory.  The  stable  aspect  of 
the  coda,  coda  decay  with  time,  can  be  interpreted  several  dif¬ 
ferent  ways  which  may  or  may  not  include  anelasticity .  This 
result  alone  has  important  implications  in  the  interpretation  of 
coda  at  shorter  periods.  It  suggests  that  homogeneous  coda  data 
sets  cannot  be  used  to  discriminate  between  coda  mechanisms  and, 
based  on  the  assumptions  inherent  in  the  models  considered  here, 
that  the  gross  distribution  of  scatterers  is  as  important  to  coda 
generation  as  is  the  density  of  scatterers. 

Ravleigh  Wave  Data 

Figure  1  shows  the  location  of  NTS  and  the  station  distribu¬ 
tion  used  in  this  study.  Long-period  digital  Rayleigh  wave  data 
were  obtained  from  the  Center  for  Seismic  Studies,  Arlington, 
Virginia,  for  a  series  of  underground  nuclear  explosions  at  NTS. 
The  data  were  previously  digitized  at  a  sampling  rate  of  1 
sample/sec  and  were  available  from  the  Center's  waveform  dar _ 
base.  Both  WWSSN  and  Canadian  network  stations  were  used.  Table 

1  lists  event  parameters  for  the  underground  explosions  and  Table 

2  lists  station  parameters  with  available  events  recorded  by  each 
station.  Station  data  were  chosen  on  the  basis  of  good  signal - 
to-noise  characteristics  and  location  within  western  North  Amer¬ 
ica  . 

Figure  2  shows  the  long-period  vertical  component  for  the 
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9/06/79  explosion  recorded  at  COR.  The  raw  data  show  a  long 
duration  coda  arriving  after  the  Airy  phase,  which  is  the  largest 
arrival  on  the  seismogram.  Narrow  band-pass  filters  of  the  data 
at  center  frequencies  of  0.05,  0.1  and  0.2  hz  show  the  develop¬ 
ment  of  coda.  Three-pole,  zero-phase,  recursive  Butterworth  fil¬ 
ters  were  used  in  all  analyses  of  this  paper.  Corner  frequencies 
for  each  of  the  respective  bandpasses  were  0.04/0.0667, 

0.08/0.12,  and  0.167/0.25  hz . 

Note  that  coda  is  best  developed  for  the  0.2  hz  bandpass  in 
Figure  2  and  that  there  seems  to  be  significant  coda  even  for 
0.05  hz  Rayleigh  waves.  Figure  3  shows  the  LPZ  component  at  LUB 
with  similar  bandpass  filtered  seismograms.  In  this  case,  noise 
dominates  the  coda  for  the  0.05  hz  bandpass.  It  is  also  seen 
that  the  0 . 2  hz  bandpass  contains  much  more  energy  before  the 
peak  of  the  envelope,  compared  to  COR,  which  is  likely  to  be  from 
body  waves  and  higher  mode  surface  waves.  This  kind  of  contami¬ 
nation  will  affect  the  interpretation  of  coda  and  is  discussed 
below.  It  was  found  that  most  stations  had  high  levels  of  noise 
at  0.05  hz  compared  to  coda  after  the  vertical  component  Rayleigh 
arrival.  Thus,  0,1  and  0.2  hz  waves  were  principally  analyzed. 

Scattering  Theories  and  Data  Analysis 

Coda  amplitude  is  modeled  in  the  time  domain  by  fitting  theo¬ 
retical  curves  to  the  bandpass  filtered  data  traces.  Assuming 
that  the  coda  is  comprised  of  scattered  Rayleigh  waves,  the 
single  scattering  model  of  Aki(1969)  is 


where,  Ac(o),t)  is  the  envelope  of  coda  amplitude,  Qc  the  coda  Q 
parameter,  u  is  circular  frequency,  and  t  the  time  from  event 
origin  time,  C  is  a  constant  which  contains  the  source  spectrum 
and  density  of  scatterers. 

An  energy  flux  model  for  Rayleigh  wave  scattering  can  be  con¬ 
structed  following  Frankel  and  Wennerberg( 1987 ) .  Assuming  radia¬ 
tion  from  a  point  source  in  a  scattering  half space,  Rayleigh 
wave  energy  will  propagate  outward  with  the  group  velocity  in  an 
expanding  cylindrical  wavefront  of  radius  r.  The  total  energy  is 
given  by 


ET  ' 


E0  +  Ec 


(2) 


where  ED  the  direct  Rayleigh  wave  energy  and  Ec  coda 
wave  energy.  Following  Frankel  and  Wennerberg' s(1987)  develop¬ 
ment  for  a  cylindrical  wave,  coda  amplitude  is  given  by 


Ack 


>,U 


V"EjT 


-cat/Q  A  -wl/Z'v't 


(3) 


where  d  is  a  scale  factor  relating  Rayleigh  wave  energy  den¬ 
sity  to  amplitude.  U  is  the  group  velocity,  Qs  is  the  Q  of  the 
scattering  medium,  Qj  is  the  anelastic  or  intrinsic  Q  of  the 
medium,  and  &  is  the  average  depth  of  Ravleigh  penetration  at 


the  cylindrical  wavefront. 

The  total  energy,  Ey .  is  found  by  correcting  the  amplitude  of 
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the  direct  Rayleigh  wave  for  scattering  and  anelastic  attenuation 
suffered  in  transit.  If  Ep'  is  the  observed  direct  wave  energy, 
then 

C*  C*  tji  Jy/Q  Ct)L  ,/Qt 

Ej  =  Eq  e  d/ws  e  d  I  .  (4) 


The  direct  Rayleigh  wave  arrival  time  is  t^-r/U.  Furthermore 


U  ]D 


2rrr  6z 


(5) 


where , 


V  (t)  dt 


't, 


(6) 


and  v(t)  is  particle  velocity.  In  practice,  a  choice  of  t^ 
and  t2  must  be  made  to  estimate  where  the  main  Rayleigh  pulse 
occurs . 

Using  relations  (4)  and  (5)  to  obtain  Ej  and  substituting 
back  into  (3)  gives 

_ _  w  i 

A  (oj,t)  =  — — e  2  IQ.  ^  -  gWl/Og) i  g-wt/2Qj  (7) 

c  i  * 


The  principal  assumptions  in  this  model  are  that  scattered  coda 
waves  are  homogeneously  distributed  behind  the  Rayleigh  wave 
front  and  that  there  are  no  conversions  from  Rayleigh  to  other 
wave  types.  The  functional  form  of  (7)  allows,  in  principle,  the 
separation  of  Qs  and  Qj .  Frankel  and  Wennerberg( 1987 )  show  that 
coda  time  decay  is  controlled  mainly  by  Qj  but  that  coda  levels 
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are  weakly  affected  by  Qs. 

I~..gstonvi33G)  studied  coda  contained  within  ideseiamic 
receiver  functions  for  Pasadena  station  and  concluded  that  yet 
another  mechanism  could  affect  the  level  and  decay  of  coda  in 
addition  to  scattering  and  anelastic  attenuation.  High  levels 
and  long  duration  of  P  wave  coda  implied  that  considerable  scat¬ 
tering  occurred  in  the  crust  and  lithosphere  under  the  station. 
Analysis  of  plane  layered  but  random  velocity  structures  situated 
over  a  homogeneous  halfspace  showed  that  theoretical  P  coda  was 
strongly  affected  by  this  gross  distribution  of  scatterers  as 
well  as  their  density.  Coda  was  seen  to  decay  quickly  with  time 
for  models  with  low  velocity  variance  compared  to  models  with 
high  velocity  variance.  Thus,  a  model  with  a  large  amount  of 
wave  scattering  would  give  an  implausible  higher  Q  according  to 
the  single  scattering  approximation.  The  physical  mechanism 
which  controlled  coda  decay  was  body  wave  radiation  into  the 
half space.  A  simple  diffusion  law  for  energy  propagation  out  of 
the  scattering  layer  was  sufficient  to  parameterize  the  effect 
for  both  synthetics  and  data. 

Numerous  theoretical  studies  (e.g.,  Munasinghe  and  Farnell , 
1973;  Martel  et  al  ,  1977;  Vidale  and  Helmberger,  1988;  McLav  j'.J-in 
and  Jih,  1988)  and  several  observational  studies  (Key  1967,1968; 
Greenfield,  1971;  Aki  and  Chouet,  1975;  Langston,  1988)  have 
shown  that  heterogeneous  structure  gives  rise  to  significant  body 
wave-surface  wave  conversion  and  interaction.  Intuitively,  it  is 
easy  to  envision  the  situation  where  significant  parts  of  Rayl¬ 
eigh  wave  energy  are  converted  to  body  waves  which  then  radiate 
away  from  the  surface  to  other  distances.  Indeed,  recent  calcu- 
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lations  (McLaughlin  and  Jih,  1988)  show  that  sv'-h  a  mechanism  is 
very  important  in  the  attenuation  of  Ihz  Rayleigh  waves.  The 
converted  body  waves  which  radiate  away  from  the  free  surface 
will  take  energy  from  the  primary  Rayleigh  wave  as  well  as  from 
the  coda  arriving  afterwards.  This  radiation  effect  could  mimic 
anelasticity  as  seen  in  Langston( 1988 ) . 

Following  Langston(1988)  an  energy  flux  model  which  incorpo¬ 
rates  radiative  diffusion  can  be  constructed  for  a  scattering 
Rayleigh  wave.  Instead  of  equation  (2),  the  total  energy  flux 
for  such  a  case  can  be  written  as 


Ec  +  ER 


(8) 


where  the  new  energy  term  Eg  represents  scattered  energy  from  the 
coda  and  the  direct  Rayleigh  wave  which  diffuses  into  the  half- 
space  from  the  layer  of  scatterers .  Since  the  total  energy  of 
the  system  is  constant, 

dEp  dE ,  dED 

— -  + — ^  + - -  —  u  .  (9) 

dt  dt  dt 


Assuming  that  the  rate  of  energy  which  leaks  out  of  a  layer 
6z  in  thickness  is  proportional  to  amount  of  coda  and  direct  wave 
energy  within  the  layer  gives 


% 

dt 


=  Y  (Ec  4  E0) 


(10) 


such  that  at  t-0,  Eg-0 .  The  proportionality  constant  is  Y  , 
the  energy  diffusivity. 
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scattering  attenuation,  equation  (12),  the  radiative  diffusion 
effect  represents  one  of  the  mechanisms  inherent  in  producing 
attenuation  due  to  scattering.  The  other  mechanism  involves 
scattering  of  direct  wave  energy  into  the  coda  field  before  coda 
energy  leaks  into  the  halfspace.  Clearly,  equation  (14)  only 

_i 

makes  physical  sense  when  0  ^  QD»  Otherwise,  the  radiative  dif- 

s  K 

fusion  mechanism  would  take  more  energy  out  of  the  direct  wave 
than  is  available  and  would  violate  the  conservation  of  energy 
equation  ( 8 ) . 

The  observed  data  at  each  station  were  normalized  following 
(6).  Figure  4  shows  the  computational  sequence  for  the  9/27/78 
event  at  ALQ .  First  the  vertical  data  are  bandpass  filtered  to 
obtain  a  trace  proportional  to  velocity.  The  trace  is  then 
squared  and  integrated.  Time  t^  was  always  taken  as  the  origin 
of  the  trace  and  t2  chosen  as  that  representing  the  main  envelope 
of  the  filtered  Rayleigh  pulse.  As  seen  in  Figure  4,  the  error 
in  this  estimate  is  of  the  order  of  10%  since  the  squared  coda 
amplitude  formed  a  minor  part  of  the  signal  in  all  cases.  The 
square  root  of  Ip  was  then  divided  into  the  trace  point  by  point. 

Individual  coda  envelopes  were  then  computed  from  the  normal¬ 
ized,  filtered  traces  by  calculating  the  amplitude  of  the  ana¬ 
lytic  signal  (Farnbach,  1975).  The  envelopes  were  then  averaged 
together  at  each  station  to  smooth  fluctuations  in  the  coda. 

The  logarithm  of  equations  (1).  (7)  and  (14)  were  then  fit  to 
the  logarithm  of  the  averaged  envelope  data  using  a  least  -  squares 
algorithm.  Analytic  partial  derivatives  of  log(C),  Qc .  Qs 
Qj  ^ ,  and  Qp ' ^  were  determined  and  used  to  construct  the  normal 
equations.  The  inversion  then  iterativelv  solved  for  these  atte- 
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nuation  coefficients  except  in  the  case  of  equation  (1)  which  was 
linear  and  with  Qc and  log(C)  obtainable  in  one  iteration. 

Inversions  using  equation  (7)  were  insensitive  to  starting 
model  and  uniformly  converged  in  all  cases  within  three  iter¬ 
ations.  However,  inversions  for  Qp'^  and  Qs"^  using  equation 
(14)  were  sensitive  to  starting  model.  If  the  starting  model  was 
too  far  from  the  true  model  (based  on  synthetic  tests),  then  par¬ 
ameter  perturbations  occasionaly  caused  Qp to  be  greater  than 
Qs'l  in  violation  of  energy  conservation.  A  strategy  was  used  in 
which  Qp*^  was  estimated  from  the  slope  of  the  log  coda  envelope 
and  then  a  small  difference  in  the  two  attenuations  was  assumed 
for  a  joint  inversion.  Convergence  in  three  iterations  for  each 
case  was  again  obtained  using  this  strategy.  Reasonably  close 
but  different  starting  models  had  little  effect  on  convergence. 

The  variances  of  the  attenuation  coefficients  were  found  by 
multiplying  the  variance  of  the  data  fit  with  the  covariance  mat¬ 
rix  of  the  least-squares  system  of  equations  (eg.,  see  Aki  and 
Richards,  1980,  page  688). 

The  amount  of  coda  fit  with  these  equations  varied  according 
to  the  data  length  of  the  available  seismograms  at  each_ station 
and  the  noise  level.  It  was  found  that  the  signal - to-noise  level 
was  poor  for  20  sec  (0.05  hz)  coda  waves  where  the  noise  was  def¬ 
ined  as  the  trace  before  the  P  wave  arrival  (see  Figure  3). 
Development  of  coda  for  0.2  and  0.1  hz  Rayieigh  waves  appeared  to 
be  stable  approximately  100  sec  after  the  Rayleigh  arrival.  The 
maximum  lapse  time  was  always  defined  by  the  shortest  signal 
length  available  from  the  digital  waveform  data  base. 

Figure  5  displays  the  coda  fits  for  0.1  hz  Rayleigh  coda. 
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Figure  6  shows  fits  for  0.2  hz  coda.  Since  each  endc  theory  fit 
the  data  equally  well  only  the  fit  for  the  s.ngle  scattering 
model  is  shown.  Figures  7  through  11  show  the  inferred  Q  parame¬ 
ters  for  each  model  and  frequency. 

Discussion 

Codas  for  both  0  1  and  0.2  hz  Rayleigh  waves  are  seen  to 
decay  linearly  on  the  log-amplitude  plotr  of  Figure  5  and  6. 

Thus,  these  waves  are  similar  in  character  to  coda  seen  in  local 
seismograms  (Aki,  1969;  1980a;  Aki  and  Chouet,  1975;  Rautian  and 
Khalturin,  1978;  Rautian  et  al  1978).  The  0.2  hz  coda  likely 
contains  scattered  energy  of  other  seismic  phases  such  as  the  S 
wave  and  higher  mode  Rayleigh  waves.  This  is  evident  in  Figure  6 
where  large  pre-Rayleigh  arrivals  can  be  seen  in  data  from  the 
more  distant  stations  before  the  peak  of  the  amplitude  envelope. 
This  contamination  is  bothersome  and  will  affect  the  Q  determina¬ 
tions  in  an  unknown  way.  However,  the  linear  coda  decay  will 
simply  be  treated  as  an  empirical  observation.  The  non-unique 
interpretation  of  its  meaning  will  be  illustrated  by  the  three 
model  fits. 

The  single  scattering  model  (Figure  7)  yields  Qc  for  both 
frequencies  that  are  comparable.  Note  the  stability  of  the 
determination  for  most  stations  falling  between  100  and  250. 

This  simple  result  quantifies  the  observation  that  the  log-coda 
slope  appears  to  be  approximately  the  same  for  any  station  in 
western  North  America.  This  observation  is  consistent  with  the 
interpretation  that  coda  consists  of  waves  scattered  over  the 
western  half  of  the  continent.  Nevertheless  there  may  indica- 
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Substitution  of  (10)  into  (9)  and  solving  the  differential 
equation  for  Ep+Ec  gives 


since  at  t-0  there  is  no  coda  energy  flux,  Ec  and  the  direct  wave 
contains  all  of  the  energy  of  the  system.  If  the  effect  of  scat¬ 
tering  attenuation  is  defined  in  the  conventional  way  (e.g., 
Frankel  and  Wennerberg,  1987)  as 


Eq  -  Ej  e 


-wt/Q, 


and  (12)  is  substituted  into  (11),  the  coda  energy  flux  is  given 
by 


Er  -  Et  L 1 


2  wt{  Qs  -  Qr)  j  e-«t/0p 


where  the  "radiative  Q" ,  Qr  is  defined  by 

r  =  w/qr  . 

Following  the  energy  flux  development  above  with  this  new 
term  and  ignoring  anelasticity  gives  the  coda  amplitude 


Ac(w,t)  = 


^  eWld/20s  [  1  •  e  Wl(  Qs  '  QR]  ]  * 


e-wt/2QR  . 


Comparing  (14)  with  (7),  it  can  be  noted  that  the  radiative 
diffusion  of  converted  Rayleigh-to-body  waves  in  the  coda  mimics 
the  effect  of  anelasticity.  Because  of  the  definition  used  for 
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cions  that  receiver  effects  are  present  particularly  for  ALQ  and 
SES  at  0.1  hz . 

Mitchell  (1975)  studied  distance -ampl itude  relationships  for 
fundamental  mode  Rayleigh  waves  from  two  explosions  in  Colorado. 
His  attenuation  coefficients  were  greater  for  western  North  Amer¬ 
ica  compared  to  the  east.  If  an  average  group  velocity  of  2.88 
km/sec  is  taken  to  represent  both  0.1  and  0.2  hz  waves  and  the 
upper  bound  of  Mitchell’s  estimates  for  attenuation  coefficients 
is  used,  values  of  about  500  and  135  are  obtained  for  anelastic  Q 
measured  at  0.1  and  0.2  hz ,  respectively.  However,  Hwang  and 
Mitchell ( 1987 )  discuss  Q  models  for  western  North  America  based 
on  more  recent  measurements  made  by  Chen(1985).  Chen's  attenua¬ 
tion  coefficients  for  0.1  and  0 . 2  hz  Rayleigh  waves  yield  average 
Q's  of  about  140  for  both  wave  frequencies.  These  empirical  val¬ 
ues  agree  quite  well  with  Qc  (Figure  7)  determined  from  Rayleigh 
coda.  This  suggests  that  the  same  attenuation  processes  act  on 
direct  Rayleigh  waves  as  well  as  on  scattered  Rayleigh  waves  in 
the  coda.  A  case  can  be  made  that  anelastic  attenuation  dom¬ 
inates  both  coda  and  direct  waves,  if  it  can  be  determined  that 
contributions  from  scattering  effects  are  small . 

Separation  of  Qs  and  Qj  is  possible  with  the  energy  flux 
model,  equation  (7).  At  0 . 1  hz  (Figure  8),  Qs  is  seen  co  be  one 
to  two  orders  of  magnitude  greater  than  Qj .  Values  for  Qj  are 
comparable  to  Qc  (Figure  7)  which  demonstrates  the  sensitivity  of 
Qj  to  the  slope  of  the  log-coda  curve.  Qs  is  primarily  obtained 
from  the  zero  lapse  time  intercept  in  the  inversion.  Qs  will  be 
sensitive  to  the  way  the  direct  Ravleigh  wave  is  normalized.  Any 
propagation  effects  which  affect  the  direct  wave  can  bias  the  Qs 
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estimate.  Focussing  or  defocussing  of  the  direct  wave  by  multi  - 
pathing  is  a  plausible  and  likely  mechanism  for  variation  in  this 
parameter  with  position  in  western  North  America  (Frankel  and 
Wennerberg,  1987).  Indeed,  the  marked  low  value  of  Qg  at  PNT  is 
consistent  with  observations  made  be  Mitchell  (1975)  on  low 
amplitude  Rayleigh  waves  observed  in  the  Canadian  Cordillera. 

The  difference  between  Qs  and  Qj  at  0.2  hz  is  somewhat  less 
(Figure  9)  than  that  at  0.1  hz  but  is  still  about  one  order  of 
magnitude  Lower  values  of  Qs  at  0.2  hz  may  imply  larger  scat¬ 
tering  effects  for  higher  frequency  Rayleigh  waves.  These 
inferred  values  may  also  imply  larger  contributions  of  other 
scattered  waves  in  the  Rayleigh  coda.  Qs  will  tend  to  decrease 
if  the  likely  contamination  of  body  and  higher  mode  surface  waves 
becomes  significant.  Note  that  Qj  for  0.1  and  0.2  hz  waves  are 
again  similar. 

Interpretation  of  the  same  data  set  with  the  radiative  diffu¬ 
sion  model  points  out  the  non-uniqueness  of  coda  interpretation 
and  the  importance  of  scattering  mechanism  and  scatterer  distri¬ 
bution  in  creating  apparent  anelastic  effects.  There  is  no 
effect  of  anelasticity  parameterized  within  the  model .  Consider¬ 
ing  the  functional  similarity  of  equations  (7)  and  (14),  it  is 
not  surprising  that  the  "radiative"  Q,  <^ ,  is  nearly  the  same  as 
Qj  in  all  cases  (Figures  10  and  11)  for  each  station.  Both 
and  Qs  are  measures  of  the  scattering  efficiency  of  the  medium 
but  are  influenced  profoundly  by  the  scatterer  geometry.  It  is 
Interesting  to  note  that  in  both  cases  (Figures  10  and  11)  both  Q's 
are  nearly  the  same.  The  interpretation  could  be  made  that  the 
rate  of  energy  scattered  out  of  the  Rayleigh  wave  is  about  equal 
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to  the  rate  at  which  it  gets  radiated  into  the  mantle  as  body 
waves.  Although  coda  is  produced  by  scattering,  most  scattered 
energy  flows  out  of  the  system  into  the  mantle. 

Although  most  stations  displayed  significant  noise  for  coda 
in  the  0.05  hz  bandpass,  signals  at  MSO  and  PNT  were  adequate  for 
inversion.  Table  3  shows  the  inferred  Q  values  from  each  type  of 
inversion.  Q  values  for  MSO  determined  from  all  three  theories 
were  comparable  to  those  found  for  0.1  hz  coda. 

Those  for  PNT,  however,  show  very  anomalous  behavior.  The 
direct  Rayleigh  wave  at  PNT  was  generally  lower  in  amplitude  com¬ 
pared  to  other  stations  at  similar  distance  and  had  relatively 
higher  coda  levels.  Qc ,  Qj .  and  Qr  all  yielded  very  low  values 
near  30.  Qs  for  the  energy  flux  and  radiative  diffusion  models 
were  also  very  low.  This  anomalous  behavior  can  be  explained  if 
structure  near  the  stations  or  along  the  path  serves  to  defocus 
the  direct  wave  arrival .  Defocussing  will  cause  a  small  observed 
direct  Rayleigh  wave.  Since  the  trace  is  normalized  to  the 
inferred  direct  wave,  coda  will  be  artificially  enhanced  causing 
apparent  low  Q's. 

There  have  been  numerous  studies  of  coda  waves  (Herraiz  and 
Espinosa,  1^87)  many  of  which  make  empirical  correlations  of  Qc 
with  tectonic  province  or  other  measures  of  anelasticity .  there 
can  be  no  doubt  that  coda  ultimately  comes  from  waves  scattered 
from  heterogeneity  in  the  lithosphere.  However,  it  is  clear  from 
this  study  that  the  interpretation  of  coda  mechanism  is  still  pro- 
blamatical  and  tied  to  theory  assumptions.  Observations  of  coda 
in  a  homogeneous  data  set  (e.g.,  S  coda  from  shallow  events  or 
surface  waves  from  explosions)  may  be  interpreted  in  several 
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plausible  ways.  These  results  and  those  of  Langston( 1 988 )  sug¬ 
gest  that  experiments  be  designed  to  constrain  the  gross  geometry 
of  scatterers  in  addition  to  measuring  coda  and  direct  wave  decay 
characteristics.  It  might  be  fruitful  to  base  data  interpreta¬ 
tions  on  accurate  numerical  models  of  elastic  wave  scattering 
where  waves  of  differing  type  sample  the  same  structure  (e.g., 
Frankel  and  Clayton,  1986).  For  example,  the  combined  observa¬ 
tions  of  teleseismic  body  waves,  local  earthquake  S  waves,  and 
regional  surface  waves  in  a  particular  region  may  be  incorporated 
into  a  self-consistent  structure  and  attenuation  model.  These 
speculations  aside,  there  is  the  need  to  accurately  address  the 
wave  scattering  mechanisms  which  come  to  play  in  shaping  the 
seismic  coda. 


Conclusions 

Coda  for  0.1  and  0.2  hz  Rayleigh  waves  propagating  from  NTS 
to  western  North  American  seismic  stations  display  exponential 
decay  with  time  in  a  similar  fashion  to  short-period  S  wave  coda 
from  local  earthquakes .  The  exponential  decay  is  similar  for 
most  stations  and  implies  coda  Q's  of  100-250  for  both  frequen¬ 
cies.  This  behavior  is  consistent  with  coda  waves  being  produced 
from  heterogeneity  distributed  across  the  western  half  of  North 
America . 

Inversion  of  stacked  coda  envelopes  using  three  fundamentally 
different  theories  yields  equally  good  fits  to  the  data.  If  the 
single  scattering  model  is  assumed,  coda  Q  is  consistent  with 
previously  determined  anelastic  Q  for  the  region.  Separation  of 
scattering  Q  and  anelastic  Q  using  Frankel  and  Wennerberg ' s ( 1987 ) 
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energy  flux  formalism  again  produces  consistent  anelastic  Q  and  a 
scattering  Q  one  to  two  orders  of  magnitude  larger  than  anelastic 
Q.  A  third  theory  which  incorporates  the  radiation  of  scattered 
surface -to- body  wave  conversions  into  a  halfspace  underlying  a 
scattering  layer  can  also  adequately  explain  coda  level  and  decay 
without  recourse  to  anelasticity . 

Although  coda  and  direct  Rayleigh  wave  behavior  is  consistent 
with  the  presence  of  significant  anelasticity  in  the  crust  of 
western  North  America,  there  is  no  diagnostic  test  that  can  be 
made  with  existing  theory  on  this  homogeneous  data  set  to  discri¬ 
minate  between  scattering  and  anelastic  effects.  The  radiative 
diffusion  model  constructed  on  the  basis  of  teleseismic  P  wave 
receiver  functions  and  modified  for  the  Rayleigh  wave  data  sug¬ 
gests  that  the  overall  geometry  of  scatterers  within  the  litho¬ 
sphere  may  be  as  important  to  the  behavior  of  the  coda  as  the  scat- 
terer  density. 
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Table  1  -  NEIS  Explosion  Parameters 


Index  Event _ 

Date 

Latitude 

!  Longitude 

mb 

I 

Sandreef 

11/09/77 

22:00:00.1 

37.07N 

116. 05W 

5.7 

2 

Panir 

8/31/78 

14:00:00.2 

37.27 

116.35 

5.6 

3 

Rummy 

9/27/78 

17:20:00.0 

37.07 

116.01 

5.7 

4 

Pepato 

6/11/79 

14:00:00.0 

37.29 

116.45 

5.5 

5 

Hearts 

9/06/79 

15:00:00.1 

37.08 

116.05 

5.8 

6 

Sheepshead 

9/26/79 

15:00:00.1 

37.22 

116.36 

5.6 

7 

Kash 

6/12/80 

17:15:00.1 

37.28 

116.45 

5.6 

8 

Tafi 

7/25/80 

19:05:00.1 

37.25 

116.47 

5.5 
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ALQ 

930 

104 

1  , 

3. 

4, 

6 

COR 

980 

326 

1 . 

5 , 

6, 

7  , 

8 

GOL 

990 

72 

1, 

3, 

4, 

5, 

6 

JCT 

1710 

112 

5, 

6. 

7 

LUB 

1380 

104 

1 , 

5, 

6, 

7  , 

8 

MSO 

1060 

10 

1 , 

3, 

4, 

5, 

6, 

EDM 

1770 

7 

1, 

2, 

3, 

4, 

5, 

PNT 

1340 

350 

all 

SES 

1500 

15 

1. 

2, 

3, 

8 
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Table  3  -  inversion  Results  for  0.05  hz  coda 


Energy  Flux  Radiative 


0„ 

Or 

0. 

Op 

Jl 

±L 

MSO 

222 

2577 

232 

213 

232 

(286/182)*  (2878/2334) 

(301/189) 

(272/175) 

(300/189) 

PNT 

29 

127 

31 

25 

31 

(30/28) 

(134/121) 

(32/30) 

(26/24) 

(32/30) 

*  Numbers  in  parentheses  denote  high  and  low  values  of  Q  based 
on  one  standard  deviation  of  error  in  the  inverted  attenu¬ 
ation. 
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■  Figure  Captions 

Figure  1:  Map  of  western  North  America  showing  VWSSN  and  Canadian 
network  stations  used  for  the  study  of  Rayleigh  coda  from 
nine  underground  nuclear  explosions  at  NTS.  Also  shown  are 
the  great  circle  paths  between  NTS  and  the  stations. 

Figure  2:  Long-period  vertical  component  at  COR  (Corvallis,  Ore¬ 
gon)  showing  the  fundamental  mode  Rayleigh  wave  and  coda  for 
the  9/06/79  explosion.  Below  the  top  data  trace  are  three 
different  bandpass  filtered  traces  showing  the  development  of 
coda  at  frequencies  of  0.05,  0.1  and  0.2  hz  Amplitudes  for 
the  filtered  traces  have  been  normalized  according  to  the 
procedure  discussed  in  the  text. 

Figure  3:  Long-period  vertical  component  at  LUB  (Lubbock,  Texas) 
for  the  9/26/79  explosion.  Same  scheme  as  Figure  2.  Note 
the  large  amount  of  energy  occurring  before  the  peak  of  the 
0 . 2  hz  bandpass.  These  arrivals  are  associated  with  the 
shear  wave  and  higher  mode  surface  waves  and  likely  contrib¬ 
ute  to  the  coda. 

Figure  4:  Example  of  the  computational  steps  involved  in  normal¬ 
izing  each  data  trace.  The  top  trace  shows  the  long-period 
vertical  component  recorded  at  ALQ  (Albuquerque.  N'M)  for  the 
9/27/78  explosion.  The  data  are  bandpass  filtered,  squared, 
and  then  integrated.  Ip  is  chosen  from  the  integrated  trace 
based  on  the  duration  of  the  fundamental  mode  wave  packet 
seen  in  the  squared  trace.  The  arrow  shows  the  pick  of  t2 
for  equation  (6).  The  error  in  normalization  is  less  than 
10%  for  all  cases  since  the  coda  is  much  smaller  than  the 


peak  amplitude  of  the  fundamental  mode  Rayleigh  wave. 
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Figure  5:  Stacked  coda  envelopes  (0.1  hz)  for  each  station  with 
the  fit  for  the  single  scattering  approximation  model.  Coda 
envelopes  were  shifted  in  amplitude  in  each  graph  by  a  factor 
of  10  for  plotting  purposes.  Stations  have  been  grouped  by 
distance  (see  '’"able  2).  The  peak  in  each  envelope  corre¬ 
sponds  to  the  direct  Rayleigh  arrival.  Note  that  the  enve¬ 
lope  slopes  are  comparable  for  all  stations. 

Figure  6:  Stacked  coda  envelopes  (0.2  hz)  with  coda  fits.  Same 
scheme  as  Figure  5. 

Figure  7:  Inferred  Qc  values  for  the  single  scattering  model. 
Error  bars  show  one  standard  deviation  in  the  parameter  esti¬ 
mate  from  the  inversion  algorithm  Triangles  are  for  0.1  hz 
coda.  Circles  are  for  0.2  hz  coda.  Note  the  clustering  of  Qc 
values  between  100  and  200. 

Figure  8:  Qj  and  Qs  determinations  for  0 . 1  hz  coda  using  the 
energy  flux  model,  equation  (7).  Note  the  change  in  scale 
for  Qs .  Scattering  Q's  are  one  to  two  orders  of  magnitude 
greater  than  anelastic  Q's.  Anelastic  Q  clusters  around  100. 

Figure  9:  Qj  and  Qs  determinations  for  0.2  hz  coda  for  the  energy 
flux  model.  Results  are  similar  to  those  at  0.1  hz . 

Figure  10:  Scattering  0,  Qs ,  and  "radiative"  Q,  ,  for  0.1  h~ 
coda  inferred  using  equation  (14).  In  this  case,  radiative 
diffusion  of  Rayleigh- to-body  wave  conversions  out  of  the 
crustal  waveguide  mimics  the  effect  of  anelasticity . 

Figure  11:  Radiative  diffusion  model  results  for  0.2  hz  coda. 
Same  scheme  as  Figure  10.  Note  that  Qs  and  are  nearlv  the 


same  for  all  cases. 
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ABSTRACT 


Pnl  waveforms  from  two  moderate-sized  eanhquakes  in  Zambia  are  used  10 
determine  an  upper  mantle  P-wave  velocity  model  for  southern  Africa.  The  events 
are:  5/15/68,  mb=  5.7,  depth=28  km;  and  12/2/68.  mh~ 5.9.  depth=6  km.  Focal  param¬ 
eters  for  these  events  are  constrained  by  previous  workers  from  telcseismic  body 
wave  inversion.  Synthetic  seismograms  are  generated  for  various  mantle  velocity 
models  using  a  wavenumber  integration  method  until  an  acceptable  fit  to  the  data  is 
obtained.  Quality  of  fit  is  measured  primarily  by  the  Pr  'PL  amplitude  ratio.  Source- 
station  geometry  also  allows  for  the  independent  sampling  of  the  upper  mantle 
beneath  the  Kapvaal-Rhodesian  craton  and  the  mobile  belt  provinces.  Synthetics  from 
a  three-layer  crust  over  hrlfspace  mantle  model  do  not  show  prominent  precursor 
arrivals  seen  in  the  data;  these  are  interpreted  to  be  P-wavcs  turning  in  the  upper 
mantle.  The  synthetics  also  give  a  too  low  PniPL  amplitude  ratio.  Synthetic:  for 
models  with  a  mantle  P-wave  velocity  gradient  of  0. 00333/sec  fit  the  cratonic  path 
data  very  well.  Since  there  is  no  indication  of  interaction  with  a  low-velocity  zone, 
this  gives  a  minimum  lithospheric  thickness  of  120  km.  A  slightly  lower  gradient  is 
indicated  for  the  mobile  belt  regions,  with  a  minimum  lithospheric  thickness  of  140 
km.  Though  the  data  set  is  small,  there  is  no  evidence  for  a  maior  low  velocity  zone 
beneath  either  province.  Different  velocity  gradients  between  the  two  provinces 
implies  different  temperature  structure,  whicti  supports  the  hypothesis  that  a  deep, 
cool  lithospheric  root  exists  beneath  the  Kapvaal-Rhodesian  craton. 
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INTRODUCTION 


Determining  seismic  velocities  in  the  uppermost  mantle  beneath  southern  Africa 
is  important  for  constraining  thermal,  petrologic,  and  tectonic  models  of  this  region. 
Data  from  other  areas  indicate  that  upper  mantle  P-wave  velocities  and  structure 
correlate  with  surface  tectonic  provinces  and  heat  flow  anomalies.  For  example, 
beneath  stable  shield  regions,  the  upper  mantle  low-velocity  zone  (LVZ)  is  absent  or 
weak,  and  P-wave  velocity  usually  increases  with  depth  (Bn:ne  and  Dorman,  1963). 
The  upper  mantle  in  younger,  active  tectonic  regions,  like  western  North  America, 
consists  of  a  relatively  thin  "lid"  (=60  km  thick)  with  a  negative  gradient  or  constant 
velocity,  overlying  a  pronounced  LVZ  (Helmberger,  1972;  Burdick  and  Helmberger, 
1978).  These  variations  arise  from  thermal  effects  related  to  the  age  of  the  most 
recent  tectonistn  in  the  region,  as  well  as  lateral  compositional  heterogeneities  in  the 
upper  mantle.  Thus,  data  on  upper  mantle  velocity  structure  yields  information  on 
thermal  and  petrologic  properties  and  the  depth  extent  of  past  tectonic  events  (Bott, 
1982). 

Efforts  to  determine  the  upper  mantle  and  crustal  structure  in  Africa  have  con¬ 
centrated  mainly  on  the  Fast  African  rift.  By  comparison,  work  in  southern  Africa 
has  been  limited.  Early  studies,  such  as  those  done  by  Gane,  et  al.  (1956),  Hales  and 
Sacks  (1959),  and  Willmore.  et  al.  (1952)  examined  crustal  structure  in  the  Transvaal 
using  refraction  experiments,  but  such  attempts  were  generally  limned  to  determining 
average  crustal  properties. 

Studies  of  deeper  structure  have  been  performed  using  long-path  surface  wave 


dispersion.  Gumper  and  Pomeroy  (1970)  determined  phase  and  group  velocities  over 
the  length  of  the  rift  valley  and  central  southern  Africa.  Though  their  study  provides 
a  useful  crustal  model  for  our  synthetic  computations,  the  upper  mantle  velocities  are 
more  representative  of  the  rift  system,  or  at  best,  an  average  between  rifted  and  more 
stable  regions.  Bloch,  et  al.  (1969)  studied  multi-mode  surface  wave  dispersion  in 
the  same  region  of  southern  Africa  examined  in  this  study.  Although  his  method  is 
most  sensitive  to  the  shear  wave  velocity  structure,  we  can  compare  our  P-wave  velo¬ 
city  model  with  their  model. 

Green  (1978)  determined  an  upper  mantle  P-wave  velocity  model  for  eastern 
and  southern  Africa  by  inverting  P-wave  travel  times.  However,  the  upper  250  km 
of  his  model  is  more  representative  of  eastern  Africa  and  the  rift  zone,  because  of  the 
source-station  distribution  used.  Only  velocities  below  this  depth  can  be  compared 
with  those  from  this  study.  Thus,  little  work  has  been  done  to  directly  determine  the 
P-wave  velocity  structure  above  250  km  depth  in  southern  Africa. 

Several  moderate  sized  earthquakes  have  occurred  on  the  African  continent 
since  the  initiation  of  the  WWSSN  network  (Fairhead  and  Girdler,  1971).  The  objec¬ 
tive  of  mis  study  was  to  use  some  of  these  earthquakes  to  directly  determine  mantle 
P-wave  velocity  structure.  We  modeled  P ^  waves  from  two  of  these  events  located 
in  Zambia,  and  recorded  at  southern  African  stations.  Source  parameters  for  these 
events  have  been  well  constrained  by  previous  workers  (Fairhead  and  Girdler,  1971; 
Wagner.  1986).  so  we  can  utilize  the  regional  seismograms  to  obtain  earth  structure. 
Synthetics  are  calculated  for  various  upper  mantle  models  using  a  wavenumber 
integrum,.,  algorithm,  and  compared  to  data  io  determine  the  upper  mantle  P-wave 


velocity  regime.  Results  are  then  compared  with  other  models  derived  from  different 
methods.  We  also  examine  the  bearing  that  our  derived  velocity  model  has  on  the 
thermal  state  of  the  upper  mantle  in  southern  Africa. 

GEOLOGY  AND  TECTONICS  OF  SOUTHERN  AFRICA 

Modern  tectonism  in  Africa  is  confined  mainly  to  the  East  African  rift  zone.  Our 
study  region  lies  south  of  this  zone,  though  several  lines  of  evidence  indicate  that  the 
rifting  is  propagating  southwestwards.  Figure  1  is  a  geologic  map  of  the  study  area, 
showing  the  mapped  crustal  provinces  and  their  ages.  The  eastern  portion  of  the  area 
is  dominated  by  the  Precambrian  Rhodesian  and  Kapvaal  cratons,  dated  at  greater 
than  2.7  Ga,  separated  by  the  slightly  younger  Limpopo  mobile  belt.  These  are 
thought  to  have  acted  as  a  tectonically  undisturbed  unit  since  the  Am^ean.  Sui round¬ 
ing  this  assemblage  are  the  mostly  younger  platform  sedimentary  rocks  of  the  mobile 
belts  such  as  the  Damara.  Zambezi,  Irumide,  and  Gariep  belts.  Kroner  (1977)  gives  a 
thorough  description  of  these  African  crustal  provinces. 

Several  geophysical  phenomena,  generally  regarded  as  evidence  of  rift  propaga¬ 
tion  southwestwards,  are  notable  in  the  region  just  north  of  our  study  area.  O  vman 
and  Pollack  (1977)  reported  an  anomalously  high  heat  flow  in  western  Zambia.  Hear 
producticn  measurements  in  the  region  indicate  that  near  surface  radioactivity  can 
account  for  only  half  of  the  observed  anomaly.  They  suggest  that  the  excess  heat 
comes  from  the  asthenosphere.  and  that  the  lithosphere  has  been  thinned  to  less  than 
60  km.  Girdier  (1975)  described  an  extensive  negative  Bouguer  anomaly  over 
Africa.  While  most  of  the  anomaly  can  be  correlated  with  surface  expressions  of  the 


Figure  1.  Geologic  province  map  of  southern  and  centra!  Africa  taken  from  Condie 
(1982).  Events  used  in  this  study  occurred  'in  wester:.  Gambia.  Note  juxtaposition  of 
old  cratons  and  younger  mobile  belts.  Our  events  occurred  in  western  Zambia,  along 
the  Dnm-ra  Fold  Belt. 
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rift,  an  arm  running  east-west  through  northern  Zambia  and  Angola  can  not.  They 
attribute  this  to  thinned  lithosphere,  with  upw-elling  of  hot,  low  density  asthenosphere. 
de  Beer,  et  al.  (1975),  noted  an  electrical  conductivity  anomaly  in  northern  Botswana 
and  Namibia,  which  they  attribute  to  a  highly  fractured,  conductive  crust  associated 
with  incipient  rifting,  since  part  of  the  anomaly  is  continuous  with  surface  expres¬ 
sions  of  the  rift.  In  Botswana.  Reeves  and  Hutchins  (1975)  detected  a  linear  seismi¬ 
city  zone  trending  NE-SW.  All  these  observations  suggest  instabilities  in  the  litho¬ 
sphere  associated  with  rift  propagation  southwestwards. 

Because  of  the  lack  of  similar  large  -scale  phenomena,  thermal  effects  related  to 
incipient  rifting  probably  have  little  influence  on  mantle  structure  in  southern  Africa. 
However,  one  interesting  observation  in  southern  Africa,  unrelated  to  rifting,  further 
motivates  this  study.  In  addition  to  determining  the  upper  mantle  P-wave  velocity 
structure,  we  hoped  to  detect  possible  differences  in  upper  mantle  structure  beneath 
the  cratons  and  mobile  belts.  This  idea  comes  from  a  paper  by  Ballard  and  Pollack 
(1987),  who  described  a  difference  in  observed  heat  flow  between  the  cratonic  and 
mobile  belt  provinces.  They  cite  reasons  for  this  as:  1)  differences  in  crustal  heat  pro- 
ducrior.  between  the  r  vo  provinces,  and  2)  diversion  of  heat  into  the  mobile  be’.s  bv 
a  deep,  cool  lithospheric  root  beneath  the  craton.  Their  models  indicate  that  at  least 
50%  of  the  difference  may  be  attributable  to  diversion  of  heat,  with  the  root  extend¬ 
ing  from  200  km  to  400  km  depth. 

Such  a  deep,  cool  root  may  have  a  velocity  structure  different  from  that  of  the 
surrounding  hotter  sub-mobile  belt  mantle,  and  this  may  be  detectable  in  seismic 
data.  Pressure  effects  and  co.  : position  being  equal  between  the  two  regions. 


velocities  in  the  hotter  upper  mantle  should  be  lower  than  in  the  cooler  lithospheric 
root  (Bott.  1982).  For  a  homogeneous  upper  mantle,  a  temperature  gradient  of  more 
than  8  to  9  K/km  would  offset  the  pressure  effect,  causing  P-wave  velocities  to 
decrease  (Anderson  and  Sammis,  1970).  Thus,  the  velocity  gradient  determination 
can  help  constrain  the  mantle  temperature  gradient.  Fortunatelv,  the  events  and  Ma¬ 
rions  used  in  this  study  are  well  situated  to  address  tins  auestion.  A  glance  at  Figure 
2  shows  that  raypaths  from  the  two  Zambia  events  cross  either  mostly  cratonic  or 
mostly  mobile  belt  provinces,  allowing  for  independent  sampling  of  the  two  regions. 
So  our  results  could  corroborate  the  deep  cratonic  root  hypothesis. 

DATA  AND  ANALYSIS  METHOD 

The  objective  of  this  study  was  to  determine  an  upper  mantle  P-wave  velocity 
structure  for  southern  Africa.  Our  approach  was  to  forward  model  regional  Pnl 
waveforms  recorded  at  several  stations  from  two  moderate-sized  earthquakes  in  Zam¬ 
bia.  Figure  2  shows  the  locations  of  the  events,  stations,  and  raypaths.  The  effects  of 
lateral  crustal  heterogeneities  on  the  PL  waveform  were  ignored;  these  would  cause 
focusing  and  defocusing,  as  well  as  scattering,  of  PL  energy.  Forward  modeling  of 
body  waves  has  not  previously  been  applied  in  this  region.  The  advantage  is  that  this 
technique  samples  upper  mantle  properties  on  a  finer  scale  than  surface  wave 


methods. 


Figure  2.  Events  and  stations  used  in  the  study.  Shown  are  raypaths  from  each  event 
fo  stations  providing  data  tor  that  event.  Region  enclosed  by  boid  line  is  the 
KapvaaJ-Rhodesian  craton.  Note  that  raypaths  to  BUL  and  PRE  cross  only  craton; 
those  to  WIN  and  SDR  cross  only  mobile  belts.  Tick  marks  are  10°  apart. 


TABLE  I:  Event  Parameters 


Event  Date 

Snake 

Dip 

Rake 

mb  Ma  t'xl 0 2"dyne  -cm  ) 

Depth  (km 

1  5/15/68 

49" 

40"  SE 

263" 

5.7  3.32 

28 

2  12/2/68 

40  °SE 

265" 

5.9  4.38 

6 

Table  I  summarizes  the  source  parameters  for  the  two  events  used,  as  deter¬ 
mined  by  Wagner  (1986)  from  moment  tensor  inversion.  Event  1  occurred  in  western 
Zambia  on  May  15,  1968.  Inversion  gave  a  nearly  pure,  high-angle  normal  fault 
mechanism  at  28  km  depth,  with  a  large  (17%)  CLVD  component,  and  a  5  sec 
source  time  function.  A  large  earthquake  this  deep  implies  a  high  strain  rate  at  depth 
lending  support  to  the  rift  propagation  hypothesis.  Event  2  occurred  350  km 
northwest  of  Event  1.  Inversion  yielded  a  nearly  pure  normal  mechanism,  with  r. 
small  CLVD  component  and  7  cec  time  function,  at  a  depth  of  6  km.  All  four  sta¬ 
tions  (PRE,  BUL,  WIN,  SDB)  contributed  data  for  Event  1,  while  only  BUL  and 
WIN  provided  data  for  Event  2. 

A  Pnl  wave  is  defined  as  the  long-period  waveform  recorded  at  regional  dis¬ 
tances  (<1700  km  or  15°)  from  crustal  earthquakes,  though  it  has  been  observed  out 
to  25"  in  some  regions.  It  includes  the  initial  P  pulse  and  lasts  until  the  shear  wave 
arrival  time.  P ^  has  been  studied  by  Wa'lace  (1983),  Shaw  and  Orcutt  0984),  and 
Helmberger  and  Engen  (1980).  Figure  3  shows  a  sample  P ^  waveform  from  an 
earthquake  in  Zambia  (this  studyg  The  name  P nI,  as  coined  by  Helmberger  and 
Engen  (1980),  implies  that  the  waveform  is  a  combination  of  two  wave  phenomena: 
a  partially  trapped  crustal  wavetrain  (PL),  and  upper  mantle  phases  (Pn  ). 


The  Pn  portion  of  the  waveform  (not  to  be  confused  with  the  mantle  head  wave 
phase  Pn)  arrives  first,  and  consists  of  P-wave  interactions  with  the  uppermost  man¬ 
tle,  such  as  the  similarly  named  mantle  head  wave,  and  upper  mantle  turning  wave 
phases  such  as  P,  pP,  and  sP  (if  these  phases  exist).  For  a  plane  layered  earth,  these 
phases  travel  with  horizontal  phase  velocities  greater  than  or  equal  to  the  P-velocitv 
at  the  top  of  titc  mantle.  In  Figure  3,  the  first  three  impulsive  arrivals  (open  arrows i 
constitute  Pn .  Thus,  because  it  is  composed  of  P-waves  that  sample  the  uppermost 
mantle,  this  portion  of  PrJ  is  the  most  important  for  inis  study. 

PL  is  the  long-period  wavetrain  that  follows  Pr  in  Figure  3.  It  was  firs: 
identified  by  Somviile  (1930),  and  has  since  been  studied  by  numerous  workers 
(Oliver  and  Major,  1960;  Oliver,  1 964;  lieimberger,  1972;  Helmberger  and  Eager:. 
1980;  Wallace  1983;  Shaw  and  Orcutt.  1984).  It  propagates  as  partially  trapped  r- 
SV  reverberations  in  the  crust,  leaking  SV  energ'"  into  the  mantle,  a  nheno’ 'snon 
termed  "leaky  mode"  propagation,  hence  explaining  why  PL  attenuates  rapidly  with 
distance.  These  reverberattens  have  ray  parameters  greater  than  the  P-wave  slowness 
of  the  uppermost  mantle.  Panicle  motion  for  the  PL  wave  is  rrograde  eliipt’cal.  Shaw 
and  Orcutt  (!984i  examined  PL  in  great  ucn.ii  using  wavenumber  uiteurav  ,  ;eve:  I- 
inc  that  it  is  most  sensitive  to  crusts’:  thickness  and  av erase  crus;::,  veiocit  .  ar.ti  that 


structure  oeneatn  tne  crus;  a 
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Figure  3.  Sample  waveform  (vertical  component!  from  a  crust;;;  earthquake  O' 
km  depth)  in  Zambia,  recorded  at  PRE  ( R  —  1 1 00  km).  Onset  of  /'r,  PL.  and  S  is 
indicated.  The  three  open-arrowed  impulsive  a; rivals  constitute  Pr.  and  propagate  in 
tiie  uppeimost  mantle.  PL  is  the  long-r>eri<xi  wavetvain  following  P and  endin'.'.  at 
me  S  amvai. 


cases,  the  tangential  motion  was  negligible,  implying  that  little  multi-pnthhtg  energy 
wa  present. 

li  generate  synthetic  seismograms  we  used  the  wavenumber  integration  routine 
from  Bark  -r  (19X-).  Tins  method  stems  from  the  solution  in  cylindrical  coordinates 
of  the  Fourier  transformed  (Helmholtz)  wave  equation.  'Hie  frequency  domain  solu¬ 
tions  for  the  P-.  FV-.  and  SH  wave  potentials  arc 

0. r  T.r  ;wi  =  £  A  (/ 0)  j  J{(kr  \  dk 

:  6 

xjir  : o«  -  V  A  •' / Q )  |  \j /,  (:.  :k  ,co)  J,(kr  )  dk 
/  0 

yjr.b. ::<j))  =  •H/0)  |  '/j(::k. w)  Ja: 

;  (i 

in  which  the  azimuthally  dependent  terms  .-.  '10)  are  the  horizontal  source  radiation 
pattern,  J/fkr)  is  die  Bessel  Fetation  of  order  /.  k  is  the  wavenumber  (=  o/c),  co  is 
angular  frequency,  c  is  horizontal  phase  velocity,  a  d  (r.8.z)  are  the  cylindrical  coor¬ 
dinates.  Notice  that  since  the  integration  is  over  real  v  wenunttVi,  this  method  w<M 
theoretically  produce  a  complete  synthetic  seismogram.  Ir.  ">ractice.  how  er.  the 
integration  is  truncated,  so  that  the  final  response  includes  on.  •  waves  over  a 
wavenumber  band  of  interest. 

The  vertically  dependent  functions  0 ;,V/,  and  ir.  the  integrand  are  computed 
using  the  Thorr.ion-ria.sk el  propagator  matrix  method  (Thomson.  1 950:  .Tasked. 


displacement-stress  discontinuities  at  the  source  depth.  Also  included  is  the  delta 
matrix  formulation  of  Duncan  (1965)  and  Thrower  (1%5)  to  avoid  precision  errors 
when  subtracting  exponentially  growing  terms  for  the  P-SV  case.  Frequeiw- 
independent  attenuation  (<5_!)  is  included  by  assuming  complex  layer  velocities. 

The  integration  over  wavenumber  is  performed  following  Apse!  (1°7“),  where  a 
quartic  polynomial  is  fit  to  the  integrand,  and  the  integration  performed  ar.alvticallv 
over  that  polynomial.  Inclusion  of  complex  velocities  allows  for  integration  over  real 
wavenumber,  since  this  moves  the  Rayleigh  poles  off  the  real  wavenumber  axis.  The 
sampling  can  be  coarse  where  the  integrand  varies  slowly,  and  finer  where  it  varies 
rapidly. 

Finally,  the  frequency  domain  Green's  functions  for  the  principal  dislocation 
..ources  are  obtained.  These  are  then  multiplied  with  the  source  spectrum,  a  WWSSN 
long-period  instrument  response,  combined  to  yield  a  synthetic  spectrum,  and  inverse 
Fourier  transformed  into  the  time  domain. 

STRUCTURE  MODELING 

All  models  for  which  synthetics  were  calculated  had  the  same  crusiJ  parame¬ 
ters.  but  different  mantle  structure.  For  f’  ,tal  model  we  chose  model  C3  from 
Pavlin  (1981),  which  is  a  modified  version  of  the  Gumper  and  Pomeroy  (1970)  cru¬ 
stal  model.  Parameters  of  this  crust  are  shown  in  Table  II.  In  his  inversion  for 
source  parameters  of  African  earthquakes.  Wagner  i !9s6;  used  a  variation  of  this 
mode’.  His  source  crustal  thickness  was  36  km  for  the  two  events  in  this  study, 
slightly  less  than  our  ?K  km  thick  crust.  For  attenuation,  an  armtrar.  G’P  and  Qk  of 


1000  and  500,  respectively,  were  assigned  to  the  crust. 


TABLE  II:  Crustal  Model 


a  (km/sec) 

P  (km/sec) 

p  (g/'cc) 

OP 

QS 

Thi Clines'-  (km) 

5.90 

3.35 

1 

o 

1000 

500 

7.0 

6.15 

3.55 

2.80 

1000 

500 

1 1.0 

6.60 

3.72 

2.85 

1000 

500 

20.0 

The  first  model  tor  which  synthetics  were  calculated  is  called  LH01.  It  has  a 
half-space  mantle  with  a  constant  P- wave  velocity  of  8.30  km/sec,  and  a  mantle  shear 
wave  and  density  structure  taken  from  Bloch,  et  al.  ( 1 969;.  Representative  values  for 
QP  and  QS  of  500  and  200,  respectively,  were  assigned  to  the  upper  mantle.  Because 
of  the  lack  of  any  mantle  velocity  gradient,  the  Pn  portion  of  P  ^  will  contain  only 
head  wave  phases:  i.e.,  waves  which  refract  along  the  Moho.  In  all  subsequent  data 
and  synthetic  comparison  figures,  the  data  is  the  thinner  trace,  the  station  name, 
azimuth,  and  distance  are  given,  and  absolute  amplitudes  are  shown  in  microns  (data 
is  the  upper  number).  Both  synthetics  and  data  have  been  low-pass  filtered  using  a 
Butterworth  filter  and  comer  frequency  of  0.8  hz. 

Of  primary  irnpotian.e  when  comparing  the  synthetics  and  data  is  the  r  .alive 
amplitudes  of  the  Pn  and  PL  portions  in  the  waveforms,  which  we  represent  as  the 
ratio  of  the  amplitude  of  the  large  arrival  at  30  sec  (referred  to  as  PnA  }  in  the  figures 
to  the  amplitude  of  the  early  PL  oscillations.  This  ratio  is  referred  to  as  the  Pn  A  :PL 
ratio.  For  a  constant  velocity  mantle,  the  PrA  :PL  ratio  in  the  synthetics  is  a 
minimum,  since  only  head  wave  and  crustai  phases  contribute  to  PnA  .  The  PnA  PL 


ratio  will  increase  for  models  with  a  positive  mantle  gradient,  because  turning  wave- 
energy  will  contribute  to  PnA  .  Hence,  the  PnA!PL  ratio  is  a  measure  of  the  upper 
mantle  influence  on  P^.  without  recourse  to  „h.,o!ute  amplitudes. 

Figure  4  shows  the  vertical  data  and  synthetics  for  Event  1  and  model  LHOL 
All  waveforms  have  been  aligned  on  the  PnA  arrival  at  ?0  sec.  Data  from  BL'L  has 
long-period  noise,  so  PnA!PL  amplitude  ratio  may  be  inaccurate.  From  calculations 
using  generalized  ray  theory,  the  first  arrival  in  these  synthetics  is  the  mantle  head 
wave  (Pn ),  and  the  second  arrival  (PnA )  is  composed  of  similar  head  waves  from 
depth  phases  like  pPmP  and  sPmP .  Notice  that,  particularly  at  BL’L  and  PRE,  the 
PnA:PL  ratio  is  underestimated  in  the  synthetics.  At  WIN  the  misfit  is  not  as  great, 
but  still  obvious.  The  synthetic  and  data  PnAiPL  ratio  at  SDB  are  very  similar,  in 
contrast  to  the  other  data.  Based  on  the  nature  of  PnA ,  these  misfits  (except  for 
SDB)  indicate  that  ther:  is  mantle  turning  wave  energy  in  the  data,  which  is  not 
included  in  the  synthetics. 

Also  observe  that  the  synthetics  do  not  show  either  of  the  dilatational  precursor 
arrivals  seen  at  PRE  (arrowed  in  Figs.  4  anu  5).  Furthermore,  the  synthetic  PnA-PL 
time  is  too  small  (i.e.,  PnA  arrives  late),  indicating  that  overall  mantle  velocities  are 
too  low.  A  positive  velocity  gradient  could  increase  the  PnA  PL  ratio  by  including 
turning  wave  energy  in  the  synthetics,  as  well  as  producing  discrete  arrivals.  Based 
on  these  criteria,  we  believe  that  a  mode!  with  a  positive  mantle  P-wave  velocity  gra¬ 
dient  is  required. 


Figure  5  shows  the  radial  components  for  Event  I  and  model  LHOi.  PnA/PL 


Figure  4.  Vertical  dam  (thin  trace:  and  synthetics  (thick)  for  Event  1,  model  LH01 
(constant  velocity  mantle  1.  Note  long-period  noise  at  BUL.  Traces  have  been  offset 
vertically  to  facilitate  comparison,  and  aligned  on  the  PnA  phase  at  30  sec.  Distance, 
azimuth,  and  absolute  ampim.nes  (microns)  are  give.  This  format  is  followed  in  ail 
subsequent  figc  es.  Notice  ,  .  nt.V-.ic  P.tA  !PL  ratio  at  all  stations  except  STB. 

Precursor  praises  (arrowed,  a.  1’r.E  are  not  produced  in  the  synthetic. 
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Figure  5.  Radiai  data  and  synthetics  for  Event  1,  model  LIKE.  Long-penod  noise  in 
the  BL'L  and  PRE  data  preclude  meaningful  f  r  A/PL  ratio  estimates,  but  the  missing 
precursors  seen  in  the  PRE  data  indicate  that  a  positive  gradient  is  needed.  WIN 
P,  \  PL  ratio  also  indicates  that  a  positive  mantle  gradient  m. tuired.  “he  F ..  A /PL 
ratio  for  both  components  of  the  SDB  data  is  useless  due  ic  instrument  p’  oleins. 


ratios  for  BUL  and  PRE  are  unreliable  because  of  the  long-period  noise.  Therefore. 


one  can  not  safely  say  that  the  synthetic  PnA/PL  ratio  at  these  stations  is  too  small. 
This  is  clearly  the  case  at  WIN  and  SDB,  though.  Furthermore,  the  precursor  phases 
at  PRE  are  missing  in  the  synthetics.  Thus,  we  amve  at  the  same  conclusion  as  for 
the  vertical  components:  neither  the  PnAIPL  ratios  nor  the  precursor  phases  at  PRE 
can  be  modeled  with  a  constant  velocity  mantle.  This  can  only  be  accomplished  by 
incorporating  a  positive  mantle  velocity  gradient  in  the  model. 

In  subsequent  comparison  figures,  the  radial  components  will  not  be  discussed  ir. 
detail,  because  the  Pnl  data  behave  as  expected  for  P-SV  waveforms.  In  general,  Pn 
phases  are  larger  on  the  vertical  than  on  the  radial,  and  the  reverse  is  true  for  the 
early  pan  of  the  PL  wave  (composed  of  shallowly  incident  energy). 

However,  the  SDB  data  presents  an  anomalous  PnA/PL  ratio  on  the  radial  com¬ 
ponent.  Here,  PnA  is  much  larger  than  PL,  while  the  vertical  PnA  is  much  smaller 
than  PL,  a  relation  not  seen  in  the  other  data.  Because  of  the  low  signal-to-noise  -;tio 
and  lack  of  calibration  pulses  on  the  original  seismogram,  we  consider  the  SDB 
PnAtPL  ratio  unreliable.  However.  Pn  will  be  useful  later  when  comparing  synthet¬ 
ics  from  different  mantle  gradient  models. 

Shown  in  Figure  6  aie  synthetics  from  model  LH01  and  data  lor  Event  2.  Radial 
components  are  not  available  for  BUL  due  to  horizontal  instrument  malfunction.  As 
in  the  Event  1  comparisons,  we  see  that  the  synthetic  PnA  IPL  ratio  is  underes¬ 
timated  at  BUL  and  WIN.  'Die  first  arrival  at  WIN  is  much  smaller  than  that  seen  in 
the  data,  and  the  interference  between  the  first  two  arrivals  in  the  data  (arrov  ed  in 


Figure  6.  Vertical  component  data  and  synthetics  for  LH01  and  Event  2.  Note 
underestimated  P.AIPL  ratio  at  both  stations,  indicating  the  need  for  a  positive  man- 
tie  gradient.  First  arrival  at  BUI.  (arrowed)  is  not  in  the  synthetics,  and  hence  is 
thought  to  be  a  turning  wave.  Note  shoulder  at  WIN  (arrowed)  from  interference  e* 
first  two  arrivals. 


Fig.  6)  is  not  modeled.  At  BL'L.  notice  that  the  first  arrival  (compressicnal.  amnved 
in  Fig.  6)  is  missing  in  the  synthetics,  just  as  the  precursor  phases  at  PRF  for  Fvent 
were  missing,  and  we  similarly  infer  this  to  he  a  turning  P-wave.  We  again  infer  thn 
a  mantle  gradient  model  would  correct  these  misfits. 

Figure  7  shows  such  a  model,  called  SACM04  (and  a  higher  gradient  model. 
SACN105,  to  be  discussed  later i.  The  crustal  parameters  are  the  same  as  those  in 
LH01;  the  mantle  now  has  an  arbitrarily  chosen  linear  P-wave  velocity  gradient  of 
0.00333/sec,  with  a  Moho  velocity  of  8.10  km/sec.  and  a  QP  and  QS  as  be: ore.  For 
die  synthetic  computations,  this  gradient  was  approximated  by  10  km  thick  layers, 
each  with  a  velocity  increase  of  0.0333  km/sec.  The  mantle  density  and  shear  velo¬ 
city  are  the  same  as  those  in  LH01,  and  /W  is  insensitive  to  the  mantle  shear  velo¬ 
city.  Shaw  and  Orcutt  (1084),  in  a  detailed  analysis  of  PL  propagation  across  Tibet, 
showed  that  PL  is  insensitive  to  mantle  structure.  Thus,  Pn  (or  PnA  )  can  med  to 
estimate  the  gradient  magnitude.  Finally,  the  earth  flattening  transformation  used  in 
Helmberger  (1973)  was  applied  to  the  entire  model. 

The  number  of  layers  needed  in  the  wavenumber  integration  to  adequately 
model  turning  waves  was  found  from  direct  P-wave  bottoming  depths  at  •he  .tries: 
station  from  ray  parameiei -distance  curves.  Such  curves  were  dete.  mined  V  _->oth 
events,  and  the  curve  for  Event  1  is  shown  in  Figure  8.  For  any  station  distance,  the 
ray’s  s!o  mess  is  the  reciprocal  of  the  velocity  at  its  bottoming  depth.  For  example, 
at  1350  km,  the  distance  of  SDB,  the  slowness  is  0.1026  sec/km;  this  gives  a  veioci 
of  8.72  km/sec,  and  reference  to  SAC.M04  gives  a  bottoming  depth  of  160  km  for  P 
Thus,  if  5  \CM04  was  the  best  fit  model  for  SDB,  it  would  be  representative  of  the 
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Figure  7.  Models  SACM04  and  SACM05,  for  a  spherically  stratified  eanh.  SACMO-t 
has  a  mantle  P-wave  velocity  gradient  of  0.00?? 3/sec;  SACM05  has  a  grad'ent  of 
0.00666/sec,  represented  by  10  km  thick  layers.  Crustal  parameters  are  taken  from 
model  C3  of  Pavlin  (1981),  and  are  a  modification  of  those  from  Gumper  and 
Pomeroy  (1970).  Mantle  shear  velocity  structure  and  density  are  taken  from  Bloch,  ct 
al.  (1969). 


Figure  8.  Rav  parameter  versus  distance  curves  tor  the  ditect  marme  P-wave  tor 
Event  i  and  SACM04  and  SACMn5.  These  are  used  to  determine  rrv  bottoming 
depths,  as  explained  in  the  text.  Note  that  ior  a  given  distance  front  the  source,  the 
rav  parameter  for  P  is  smaller  lor  SACMU5  than  for  SmCMO-J.  Hence  I  bottoms 
deeper  in  the  higher  gram-cut  model  ior  a  given  station. 


average  structure  down  to  160  kir»  depth.  If  there  is  no  evidence  in  the  data  of 
interaction  with  an  LVZ.  this  also  constrains  the  minimum  lithospheric  thickness. 


assuming  that  an  LVZ  marks  the  base  of  the  lithosphere. 

The  vertical  synthetics  and  data  for  SAC  MO 1  and  Event  1  are  shown  in  f;.  :ure 
Recall  that  the  mantle  gradient  should  cause  the  PnA  VI.  ratio  to  increase  hy  con¬ 
tributing  turning  wave  energy  to  the  waveforms.  Notice  that  the  A .  I'L  ratio  at 
BUL  has  increased  only  slightly.  This  station  is  too  close  for  the  gradient  to  have  a 
large  effect,  and  the  data  PnA  amplitude  may  be  increased  by  low-frequency  noise. 

At  PRE.  the  first  precursor  phase  is  clearly  visible  in  the  synthetics,  its  polarity  and 
amplitude  being  equivalent  to  that  ir.  the  data.  This  confirm-;  that  the  first  precursor  in 
the  data  is  the  direct  mantle  P-wnve  because  it  is  the  first  arrival,  and  its  polarity  and 
amplitude  are  consistent  with  that  of  a  synthetic  P-wave  with  a  take-off  angie  (from 
ray  parameter-distance  curves)  in  the  dilatations!  quadrant  of  the  focal  sphere.  Based 
on  travel-times,  the  energy  primarily  responsible  for  boosting  PnA  is  sP.  Due  to  its 
small  free  surface  reflection  coefficient  for  this  slowness  (-0.1),  and  near-nodal  initial 
amplitude,  pP  is  not  seen,  and  so  we  arc  left  without  an  explanation  for  the  second 
precursor  phase.  Explanations  for  this  phase  will  be  discussed  later.  However,  we 
still  conclude  from  the  PILE  data  that  SACM04  is  a  better  estimate  of  the  mantle 
structure  beneath  the  craton  than  is  LH01. 


At  WtN  in  Figure  9,  notice  that  the  synthetic  PnA>FL  ratio  is  overestimated. 
The  synthetic  mantle  P-wave  (at  IS  sec!  at  SDB  much  larger  than  that  in  the  data, 
which  is  near  noise  level.  Together  these  imply  that  more  LUTrunti  wave  energy  is 
present  in  the  synthetics  than  is  actualh  in  the  data.  Thus,  judging  from  the  ■  erticai 


Figure  9.  Vertical  component  data  and  synthetics  for  SACM04  and  Event  1.  PRE 
synthetics  now  contain  the  first  precursor  arrival--  the  direct  mantle  P-wave.  The 
second  precursor  arrival  cannot  be  modeled.  Synthetic  PnAiPL  ratio  at  BUL  is  too 
small;  that  at  PRE  is  nearly  correct;  at  WIN.  it  is  too  large.  This  implies  that  the 
SACM04  gradient  is  too  high  to  represent  the  mooiie  bel;  upper  mantle,  whiie 
approximately  correct,  if  no  iow,  for  the  cratonic  upper  mantle. 


components,  the  gradient  in  SACM04  is  too  high,  although,  as  seen  from  the  discus¬ 
sion  for  LH01  syntheucs  for  WIN,  a  small  gradient  is  required  to  model  data  from 
the  mobile  belt  regions.  Recall  that  paths  to  WIN  and  SDB  lie  along  the  mobile  belts. 

Radial  components  for  SACM04  and  Event  1  are  given  in  Figure  10.  The  first 
precursor  phase  in  the  PRE  data  is  clearly  present  in  the  synthetic,  as  it  was  on  the 
verdcal  component,  but  the  second  precursor  phase  remains  missing.  At  WIN  the 
synthetic  PnA/PL  ratio  is  slightly  overestimated.  The  initial  P-wave  arrival  at  SDB, 
predicted  by  the  synthetics,  is  either  not  in  the  data,  or  too  near  noise  level.  These 
observations  reiterate  the  conclusions  drawn  from  the  vertical  components. 

Figure  11  displays  the  vertical  data  and  synthetics  for  model  SACM04  and 
Event  2.  Clearly,  the  fit  has  been  improved  using  the  gradient  model.  The  shoulder 
produced  by  the  interference  between  the  first  two  arrivals  at  WIN  (P  and  sP)  is 
reproduced.  We  also  notice  that  the  synthetic  PnA/PL  ratio  at  both  BUL  and  WIN 
has  increased,  though  by  too  much  at  WIN  and  not  enough  at  BUL.  This  indicates 
that  SACM04  overestimates  the  actual  gradient  sampled  by  the  WIN  data,  and 
underestimates  that  sampled  by  the  BUL  data. 

The  BUL  synthetic  predicts  a  dilatational  first  arrival  (arrowed  in  Fig.  11), 
representing  the  mantle  P-wave.  This  polarity  is  consistent  with  the  take-off  angle  of 
P  from  the  focal  sphere  in  the  model.  Therefore,  the  precursor  phase  at  BUL  cannot 
be  interpreted  as  a  direct  P-wave  in  our  structure  model.  Furthermore,  allowable 
changes  in  source  parameters  (±4°  in  strike,  dip,  and  rake)  cannot  produce  such  an 
arrival.  Examination  of  short-period  teleseismic  records  shows  no  evidence  of  a 
small  magnitude  foreshock,  ruling  out  such  an  explanation.  Speculation  on  the  nature 
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Figure  10.  Radial  component  data  and  synthetics  for  SACMC'4  and  Event  1.  Long- 
period  noise  in  the  BUL  and  PRE  data  precludes  accurate  PnAIPL  ratio  estimates. 
However,  the  first  precursor  at  PRE  is  present,  though  not  the  second.  PnA  tPL  ratio 
at  WIN  is  overestimated  by  the  synthetic  anc  me  P-wave  arrival  at  SDB  (at  18  sec) 
predicted  by  the  synthetic  is  missing,  showing  mat  this  gradient  is  too  high. 
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Fisrure  11.  Vertical  component  data  and  synthetics  for  SACMCM  and  Event  — 
PnA/PL  ratio  at  B1JL  has  increased,  though  not  enough.  The  synthetic  predicts  a 
dilatatic  nal  first  arrival;  a  compression  1  arrival  is  seen  in  the  data.  This  arrival  can¬ 
not  be  explained  by  a  mantle  gradient.  The  PnAlPL  m.io  at  'V1NT  ’■*  overestimated 
slightly,  and  the  shoulder  is  present  in  the  unfiltered  synthetic. 


of  this  arrival  will  be  diverted  until  later,  but  its  occurrence  does  riot  affect  our 
requirement  for  a  mantle  gradient. 

In  an  attempt  to  further  increase  the  PnAlPL  ratio  at  BUL  for  both  events,  and 
to  see  how  well  we  could  constrain  the  gradient  magnitude,  we  computed  synthetics 
for  model  SACM05,  which  has  a  gradient  double  the  magnitude  of  that  in  SACM04. 
SACM05  is  also  shown  in  Figure  7.  Mantle  shear  velocity,  density,  and  attenuation  is 
the  same  as  in  SACM04.  For  this  higher  gradient  model,  we  expect  the  synthetic  P 
and  PnA  amplitudes  to  increase  further  relative  to  the  PL  wave  at  all  stations. 

Figure  12  displays  the  vertical  data  and  synthetics  from  model  SACM05  and 
Event  1.  At  BUL,  the  larger  gradient  has  boosted  PnA  slightly,  so  that  the  PnA/PL 
ratio  is  close  to  the  data  ratio.  However,  at  PRE  the  PnA  phase  is  now  2.5  times  as 
large  as  the  PL  wave:  a  much  larger  ratio  than  in  the  data  (PnA  IPL  =1.66).  Further¬ 
more,  the  synthetic.  PL  wave  arrives  late  at  all  stations,  because  the  PnA-PL  time  is 
larger  than  in  the  data;  i.e.,  PnA  arrives  too  early,  implying  too  high  mantle  veloci¬ 
ties.  Thus,  PRE  indicates  that  this  gradient  is  too  large;  SACM04  synthetics  fit  all  but 
the  BUL  data  better.  The  Pn  A  amplitude  at  BUL  may  be  artificially  increased  by 
long-period  noise  (apparent  in  the  fust  20  secs),  its  true  amplitude  being  1 than 
that  shown.  Alternatively,  an  increase  in  the  gradient  at  shallow  depth  in  SAC.M04, 
where  the  energy  arriving  at  BUL  turns,  might  increase  PnA  without  requiring  a 
higher  overall  gradient,  which  is  precluded  by  the  PRE  data.  The  latter  hypothesis 
was  no:  modeled,  though  the  Evem  2  BUL  data  suggest  that  this  explanation  is 
correct. 


Event  1  Vertical  SACM05 
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Figure  12.  Vertical  component  data  and  synthetics  for  SACM05  and  Event  1. 
PnA/PL  ratio  at  WIN  is  greatly  overestimated,  and  the  large  dilatational  P- wave  at 
SDB  (18  sec)  is  not  in  the  data.  The  PnAIPL  ratio  at  BUI-  is  close  to  the  data  ratio. 
At  PRE,  the  synthetic  PnA/PL  ratio  is  slightly  too  large,  indicating  that  the  noise  at 
BUL  is  causing  an  erroneous  ratio.  Notice  that  the  PnAIPL  ratio  is  not  as  highly 
overestimated  at  PRE  as  at  WIN. 
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Obviously,  the  synthetic  PnAIPL  ratio  at  WIN  in  Figure  12  is  greatly  exag¬ 
gerated,  as  we  expected  from  the  SACM04  synthetics,  which  already  overestimated 
this  ratio.  At  SDB,  the  very  large  direct  mantle  P-wave  in  the  synthetics  is  totally 
absent  from  the  data.  Thus,  the  gradient  in  SACM05  is  definitely  too  large  for  both 
cratonic  and  mobile  belt  paths,  and  that  in  SACM04  appears  to  give  the  best  fitting 
synthetics  for  Event  1  cratonic  path  data.  Radial  components  in  Figure  13  yield  the 
same  conclusion. 

Figure  14  shows  the  vertical  component  data  and  synthetics  from  SACM05  and 
Event  2.  Clearly  the  gradient  in  this  model  is  too  high  for  the  WIN  data,  as  evi¬ 
denced  by  the  very  large  synthetic  PnAlPL  ratio  at  WIN,  as  well  as  the  extreme 
absolute  amplitude  difference.  At  BUL  though,  the  fit  is  good,  except  for  the  missing 
compressional  precursor  arrival.  There  is  no  long-period  noise  here  to  induce  ampli¬ 
tude  errors  as  there  is  for  the  Event  1  BUL  data.  The  data  and  synthetic  PnAlPL 
ratios  are  very  close,  though  for  the  Event  1  data  this  model’s  synthetics  overes¬ 
timated  the  PnAIPL  ratio  at  PRE.  For  raypaths  as  close  together  (azimuthallv)  as 
those  to  BUL  and  PRE,  SACM05  could  not  be  the  best  fitting  model  for  Event  2  but 
not  for  Event  1.  This  could  be  explained  if  the  BUL  data  could  be  fit  with  a  model 
similar  to  SACM04  but  with  a  higher  gradient  at  shallow  depth,  as  discussed  rbove. 
This  was  not  modeled,  however,  as  it  does  not  affect  the  quantitative  conclusion 
about  the  average  gradient,  which  is  best  constrained  by  the  PRE  and  WIN  data. 
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Figure  13.  Rad>'al  component  data  and  synthetics  for  SACM05  and  Even:  1.  Long- 
period  noise  at  BUL  and  PRE  precludes  measurement  of  PnA/PL  ratios.  Note  large 
P  arrival  at  SD3,  missing  in  the  data,  and  the  too  large  synthetic  PnAlPL  ratio  at 
WIN,  indicating  that  the  gradient  is  too  large. 
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Figure  14.  Vertical  component  data  and  synthetics  for  SACM05  and  Event  2.  The 
PnA/PL  ratio  at  BUL  is  very  similar  to  data  ratio,  but  this  gradient  is  too  high,  as 
shown  by  the  PRE  data,  indicating  a  steeper  gradient  than  SACM04  at  shallow  depth. 
The  first  arrival  in  the  data  is  also  missing,  indicating  more  complicated  structure 
than  a  simple  gradient.  The  PnAlPL  ratio  at  WIN  is  greatly  overestimated  in  the 
svntht  ,.c. 
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DISCUSSION 


Comparison  of  the  data  and  synthetics  above  indicates  that  there  is  definitely  a 
positive  P-wave  velocity  gradient  in  the  upper  mande  beneath  southern  Africa.  Evi¬ 
dence  for  this  is  two-fold.  First,  and  most  convincingly,  the  appearance  in  the  data  of 
mantle  turning  wave  arrivals  requires  a  positive  gradient.  This  phase  is  totally  lacking 
in  the  model  with  a  constant  mande  velocity  (LH01).  Secondly,  a  larger  PnA/PL 
ratio  in  the  data  than  in  the  LH01  synthetics  implies  some  amount  of  mande  turning 
wave  energy  in  the  seismograms. 

The  data  also  indicate  that  the  mantle  gradient  is  probably  slightly  lower 
beneath  the  mobile  belt  regions  (WEN-SDB  paths)  than  beneath  the  cratonic  regions 
(BUL-PRE  paths).  This  is  based  on  the  observation  that  the  gradient  model  SACM04 
synthetics  fit  the  PRE  data  well,  while  they  overestimate  the  PnA/PL  ratio  at  WIN 
for  both  events.  The  higher  gradient  model,  SACM05,  overestimates  the  PnA/PL 
ratio  at  all  stations  except  EUL  for  both  events,  indicating  that  this  gradient  is  too 
high.  We  therefore  deduce  that  the  gradient  is  slighdy  less  then  0.00333/sec  beneath 
the  mobile  belts,  and  approximately  0.00333/sec  beneath  the  cratons,  though  the  data 
cannot  constrain  the  gradient  this  accurately.  This  difference  conveys  important  infor¬ 
mation  about 

the  tectono-thermal  state  of  the  region,  as  discussed  below. 

The  criteria  used  to  infer  the  gradient  and  its  magnitude  are  affected  by  errors  in 
the  source  parameters.  Hence,  an  examination  of  how  amplitudes  of  P,  PnA  ,  and  PL 
are  affected  by  these  factors  is  required.  Error  bounds  on  the  moment  tensor  elements 
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are  given  by  Wagner  (1986).  Since  these  do  not  relate  directly  to  the  double  couple 
parameters,  we  determined  strike,  dip,  and  rake  for  various  moment  tensors  to  deter¬ 
mine  error  bounds  on  these  quantities.  A  conservative  estimate  for  all  parameters  is: 
±3°  for  Event  1,  and  ±4°  for  Event  2.  Of  particular  interest  is  the  sensitivity  of  the 
Event  1  PRE  and  WIN  synthetics  to  these  errors,  since  these  data  best  constrain  the 
gradient  magnitude.  We  found  that  errors  in  the  dip  had  the  largest  effect  on  the 
waveform,  but  that  this  effect  was  negligible  with  regards  to  the  PnA/PL  ratio. 

Strike  and  rake  had  similarly  small  effects  on  the  synthetics,  so  that  our  conclusions 
remain  unaffected  by  the  allowable  range  in  source  parameters. 

For  comparison  with  other  velocity  models,  and  for  delimiting  the  depths  over 
which  our  model  SACM04  applies,  we  summr  ize  mantle  P-wave  bottoming  depths 
at  the  appropriate  stations.  The  P-wave  arrival  at  PRE  bottoms  at  120  km  depth 
beneath  the  Kapvaal-Rhodesian  craton.  The  P-wave  at  WIN,  100  km  further  from  the 
source,  bottoms  at  140  km  beneath  the  Damara  Fold  Belt. 

While  they  do  not  affect  our  conclusions  about  the  average  gradient,  certain 
aspects  of  the  data  could  not  be  reproduced  in  the  synthetics.  The  most  obvious 
misfit  is  that  neither  the  second  precursor  phase  at  PRE  for  Event  1,  nor  the  first 
arrival  at  BUL  for  Event  2  could  be  reproduced  with  our  models  (see  Figs.  9  and 
11).  Both  of  these  stations  lie  on  the  craton.  Something  more  complicated  than  a 
monotonic  gradient  must  be  responsible  for  these  discrepancies,  because  allowable 
changes  in  source  parameters  cannot  produce  these  arrivals.  The  phases  could  be  due 
to  small  foreshocks,  but  examinatiOh  of  short  period  teleseismic  and  regional  records 
indicated  no  such  events.  Perhaps  reflection  from  or  refraction  through  a  high- 
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velocity  or  high-velocity  gradient  layer  beneath  the  craton  is  responsible;  testing  this 
would  require  further  modeling. 

The  P-wave  velocity  model  of  Bloch,  et  al.  (1969)  contained  no  high-velocity 
zones,  though  their  methods  were  not  sensitive  to  the  P-wave  velocity  regime,  nor  to 
thin  layers.  Other  workers  have  reported  high-velocity  layers  a:  shallow  depth  in  the 
lithosphere.  Him,  et  al  (1973)  described  relatively  high  velocity  refractors  at  55  and 
80  km  depth  in  France  from  a  long  (1000  km)  refraction  profile.  A  compilation  of 
similar  long-range  refraction  results  from  many  different  areas,  showing  evidence  of 
fine  layering  in  the  lithosphere,  is  given  by  Fuchs  (1979).  Proposed  explanations  for 
this  phenomenon  include  anisotropy  (due  to  preferred  alignment  of  olivine  crystals) 
within  the  layers,  and  lateral  petrologic  and  chemical  heterogeneities.  Bamford  (1977) 
has  demonstrated  a  7%  anisotropy  in  Pn -velocity  in  western  Germany,  though  a  sub¬ 
sequent  study  (Bamford,  et  al.,  1978)  in  northern  Britain  failed  to  reveal  any 
significant  anisotropy.  Thus,  such  layering  may  not  be  ubiquitous,  and  more  data  is 
required  to  conclusively  identify  and  study  it  in  southern  Africa. 

Precise  modeling  of  other  aspects  of  the  data  would  require  a  small  modification 
to  the  gradient  model.  Recall  that  the  SACM04  synthetics  consistently  underestimated 
the  PnA/PL  ratio  at  BUT  for  both  events  (Figs.  9  and  11).  This  was  true  even  for 
Event  2,  which  contained  no  low-frequency  noise.  Model  SACM05  produced  better 
fitting  synthetics  at  BUL,  but  not  at  PRE  for  Event  1,  and  the  latter  constrains  the 
gradient  better.  Since  BUL  is  nearer  the  source,  this  could  indicate  that  a  higher  gra¬ 
dient  is  present  at  shallow  depth,  where  the  rays  to  BUL  turn.  Such  a  modification 
would  increase  the  PnAIPL  ratio  at  BUL  without  affecting  further  stations.  Since  our 
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purpose  was  to  determine  the  average  gradient,  this  aspect  was  not  modeled. 

Synthedc  mantle  turning  wave  amplitudes  are  also  affected  by  om  choice  of 
upper  mande  QP  and  QS  values.  Higher  Q  values  would  require  a  smaller  velocity 
gradient  to  produce  similar  amplitudes,  though  relative  travel-times  would  also  be 
changed.  When  compared  with  published  Q  values  (Anderson,  et  al.,  1965;  Anderson 
and  Han,  1978;  Patton,  1980),  our  values  are  approximately  twice  as  large.  However, 
these  studies  generally  find  an  averaged  Q  over  broad  regions;  no  attenuation  studies 
have  been  undenaken  in  southern  Africa.  Therefore,  our  values  may  be  high,  but  not 
unreasonable  for  an  old,  stable  platform  and  cratonic  region.  If  a  subsequent, 
independent  Q  survey  should  yield  a  lower  Q  value,  it  could  imply  a  slightly  higher 
mande  gradient  than  determined  here. 

A  positive  mande  P-wave  velocity  gradient  has  implications  for  the  southern 
African  upper  mande  thermal  regime.  Since  our  data  reveal  that  the  gradient  is 
smaller  beneath  mobile  belts  than  beneath  cratons,  this  study  supports  the  deep  lithos¬ 
pheric  root  hypothesis.  Assuming  pressure  changes  with  depth  are  roughly  equal 
between  the  two  provinces,  and  assuming  similar  petrologies,  the  higher  velocity  gra¬ 
dient  indicates  lower  temperatures  and  temperature  gradients  in  the  cratonic  root.  The 
nominal  model  used  by  Ballard  and  Pollack  (1987)  predicts  this  difference  in  thermal 
structure.  It  indicates  that  at  a  depth  of  100  km,  the  mobile  belt  mande  is  about 
600°  C  hotter  than  the  cratonic  root,  and  the  temperature  gradient  is  higher  in  the 
mobile  belt  lithosphere  down  to  100  km.  Temperature  gradients  below  100  km  are 
higher  in  the  cratonic  root  than  in  the  adjacent  asthenospheric,  mobile  belt  mantle, 
but  between  these  two  regions  there  is  a  major  rheologic  difference.  The  temperature 
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profiles  do  not  merge  until  800  km  depth. 

A  lower  limit  to  the  depth  of  the  root  can  be  ascertained  from  the  bottoming 
depths  of  the  direct  mantle  P-waves  at  PRE  for  Event  1.  The  bottoming  depth  from 
SACM04  for  the  P-wave  at  PRE  is  120  km,  and  there  was  no  evidence  for  a  P-wave 
LVZ.  Ballard  and  Pollack  (1986)  estimate  that  the  root  may  extend  anywhere  from 
200  to  400  km  depth,  so  our  conclusion  is  well  within  the  lower  limit.  Brune  and 
Dormann  (1963),  using  surface  wave  dispersion,  did  not  detect  a  P-wave  LVZ 
beneath  the  Canadian  shield  above  400  km.  Like  the  Rhodesian-Kapvaal  craton,  the 
Canadian  shield  also  has  relatively  low  heat  flow.  In  northwestern  Eurasia.  Given  and 
Helmberger  (1980),  using  body  waves,  did  detect  a  small  LVZ  for  P-waves,  50  km 
thick,  beginning  at  150  km  depth.  Thus,  our  data  may  not  sample  deep  enough  to 
detect  an  LVZ,  but  it  does  yield  a  minimum  lithospheric  thickness  in  the  craton. 

Fairhead  and  Reeves  (1977)  have  estimated  lithospheric  thickness  over  most  of 
the  African  continent  by  combining  teleseismic  delay  times  and  Bouger  anomalies. 
Their  map  shows  a  thickness  of  150-175  km  beneath  BUL  and  PRE,  which  is  con¬ 
sistent  with  our  derived  minimum  lithospheric  root  thickness  of  120  km,  and  with  the 
Ballard  and  Pollack  (1987)  study. 

Their  map  shows  thinner  lithosphere  in  the  mobile  belt  regions,  the  thickness 
beneath  WIN  and  SDB  being  less  than  150  km.  The  direct  P-wave  at  WIN  for 
SACM04  for  Event  1  bottoms  at  140  km;  however,  since  we  believe  that  the  gradient 
is  slightly  less  than  that  in  SACM04,  the  ray  actually  bottoms  shallower.  This  depth 
is  close  to  the  lithospheric  thickness  in  the  region  given  by  Fairhead  and  Reeves 
(1977),  so  the  actual  rays  could  be  encountering  a  low  velocity  cone,  an  effect  which 
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we  did  not  model.  However,  there  was  no  hard  evidence  from  our  small  data  set  of 
interaction  with  an  LVZ,  so  the  lithosphere  is  somewhat  less  than  140  km  thick 
beneath  the  mobile  belts.  Unfortunately,  the  SDB  data,  which  would  have  given  a 
deeper  sampling  of  the  structure,  proved  unusable. 

Figure  15  compares  our  velocity  models  with  those  from  previous  studies. 
Green’s  (197S)  model  was  obtained  from  inversion  of  P-wave  travel  times  for  earth¬ 
quakes  throughout  eastern  and  southern  Africa,  the  Red  Sea,  and  Gulf  of  Aden. 
Because  of  source-station  geometry,  the  upper  250  km  of  1m  model  is  more  represen¬ 
tative  of  the  rift  zone  upper  mantle,  or  at  best  an  average  of  stable  and  active  upper 
mantle  regions.  Only  velocity  values  below  this  depth  apply  to  southern  Africa.  It 
shows  an  overall  slight  positive  gradient,  but  the  values  are  considerably  less  than 
those  in  SACM04.  This  is  probably  a  result  of  the  different  regions  examined  in  the 
two  studies,  and  our  models  do  not  extend  deep  enough  to  compare  with  Green’s 
below  250  km. 

Figure  15  also  shows  the  P-wave  velocity  model  of  Gumper  and  Pomeroy 
(1970),  obtained  from  surface  wave  dispersion.  Notice  that  it  agrees  closely  with  that 
of  Green  (1978),  probably  because  these  models  both  average  structure  between  the 
rift  zone  and  more  stable  regions. 

Finally,  we  compare  our  models  with  that  of  Bloch,  et  al.  (1969).  Their  study 
region  is  nearly  identical  to  ours,  and  there  model  does  not  average  between  major 
tectonic  provinces.  Our  model  SACM04  and  the  Bloch  model  coincide  closely  for 
most  of  the  former’s  depth  range.  Model  SACM05  diverges  rapidly  from  theirs 
below  70  km  depth,  and  we  have  shown  that  this  gradient  is  too  high.  SACM04 
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Figure  15.  Comparison  of  models  SACM04  and  SACM05  with  models  of  Gumper 
and  Pomeroy  (1970),  Green  (1978),  and  Bloch,  et  al.  (1969).  Models  are  for  a  spher¬ 
ically  stratified  earth.  The  first  two  models  are  more  representative  of  the  tectonically 
active  regions.  Note  that  SACM05  quickly  diverges  from  al!  other  models.  The 
Bloch  model  ag*-es  closely  with  out  mode:  SACM04.  resulting  from  the  two  studies 
examining  the  same  region..  SACM04  is  denned  by  the  PRE  data  to  a  depth  of  120 
km  beneath  the  craton. 
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represents  structure  at  least  down  to  120  km  depth  beneath  the  cratons.  and  probably 
deeper.  Beneath  the  mobile  belts,  the  actual  gradient  is  less  than  in  SACMO-t,  and  is 
probably  defined  to  a  depth  somewhat  less  than  140  km.  A  lower  gradient  beneath 
the  mobile  belt  regions  would  also  agTee  closely  with  Bloch's  model. 

To  determine  if  we  could  distinguish  between  a  gradient  and  a  velocity  jump,  as 
shown  in  the  Bloch,  et  al.  (1969)  model  at  70  km  depth,  we  computed  synthetics  tor 
their  model.  The  velocity  jump  model  synthetics  did  not  match  the  PRE  data  for 
Event  1,  and  so  we  conclude  that  a  gradient  is  more  appropriate  than  a  velocity  jump. 
Surface  wave  dispersion  cannot  resolve  a  velocity  gradient,  so  it  represents  it  as  a 
velocity  jump.  However,  a  model  with  a  small  velocity  jump  separating  two  layers 
having  a  positive  velocity  gradient  was  not  examined,  and  could  possibly  explain  the 
unmodeled  precursor  phase  at  PRE  and  BUL. 

CONCLUSIONS 

The  data  from  this  study  indicate  that  there  is  a  positive  P-wave  velocity  gra¬ 
dient  in  the  upper  mantle  beneath  southern  Africa.  Synthetic  seismograms  were  gen¬ 
erated  using  a  wavenumber  integration  routine,  and  compared  with  Pni  waveform 
data  from  two  Zambia  earthquakes.  Previous  workers  have  determinec  the  source 
parameters  for  these  events,  thereby  allowing  modeling  of  structure  effects.  The 
source-station  distribuuon  allows  for  examination  of  structure  beneath  cratons  and 
m-  bile  belts  independently,  to  detect  differences  in  upper  mantle  structure. 

Based  on  the  existence  in  the  data  of  mantle  turning  wave  phases  P  and  sP.  and 


102 


their  amplitudes  relative  to  PL,  the  gradient  magnitude  could  be  estimated.  A  gradient 
of  approximately  0.00333/sec  is  indicated  for  the  average  velocity  structure  beneath 
the  Kapvaal -Rhodesian  craton  (model  SACM04).  Some  of  the  data  suggest  that  the 
gradient  may  be  steeper  at  shallower  depths.  From  P-wave  bottoming  depths,  this 
gradient  extends  at  least  to  120  km,  which  is  also  a  lower  limit  on  lithospheric  thick¬ 
ness  beneath  the  craton.  This  conclusion  is  consistent  with  lithospheric  thickness  stu¬ 
dies  from  regional  gravity  and  teleseismic  P-wave  delay  times. 

Based  on  similar  criteria,  a  smaller  positive  gradient  must  exist  beneath  the 
mobile  belts.  A  magnitude  slightly  less  than  0.00333/sec,  but  still  non-zero,  is  con¬ 
sidered  likely.  For  model  SACM04  the  bottoming  depth  of  P-waves  at  WIN,  the 
furthest  mobile  belt  station  showing  no  evidence  of  a  LVZ,  is  140  km.  Lithospheric 
thickness  studies  give  a  thickness  of  less  than  150  km  for  the  mobile  belt  regions. 

Our  inability  to  model  some  precursor  arrivals  at  PRE  and  BUL  with  a  mono- 
tonic  velocity  gradient  suggests  that  there  may  be  high-velocity  layers  in  the  cratonic 
upper  mantle.  These  precursors  only  occur  at  cratonic  path  stations,  and  more  data  is 
needed  to  determine  the  nature  and  cause  of  this  phase. 

Our  conclusion  of  a  higher  velocity  gradient  beneath  the  cratonic  province  sup¬ 
ports  the  hypothesis  that  a  deep,  cool  lithospheric  root  exists  underneath  this  region. 
Our  data  also  place  a  minimum  depth  extent  on  the  root  of  120  km.  The  root  w'as 
originally  postulated  as  a  mechanism  for  heat  diversion  to  explain  lower  values  of 
observed  heat  flow  (Ballard  and  Pollack,  1987)  in  the  craton  than  in  the  mobile  belts. 
Tectonically,  this  implies  that  the  cratonic  lithosphere  has  remained  largely 
unaffected  by  post-Archean  heating  events.  Similar  data  sets  could  be  used  elsewhere 


to  detect  upper  mantle  velocities  and  heterogeneities  between  tectonically  old  and 
young  regions. 
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